int inner_solve_dprimme(double *x, double *r, double *rnorm, 
   double *evecs, double *evecsHat, double *UDU, int *ipivot, 
   double *xKinvx, double *Lprojector, double *RprojectorQ, 
   double *RprojectorX, int sizeLprojector, int sizeRprojectorQ, 
   int sizeRprojectorX, double *sol, double eval, double shift, 
   double eresTol, double aNormEstimate, double machEps, double *rwork, 
   int rworkSize, primme_params *primme) {

   int i;             /* loop variable                                      */
   int workSpaceSize; /* Size of local work array.                          */
   int numIts;        /* Number of inner iterations                         */
   int ret;           /* Return value used for error checking.              */
   int maxIterations; /* The maximum # iterations allowed. Depends on primme  */

   double *workSpace; /* Workspace needed by UDU routine */

   /* QMR parameters */

   double *g, *d, *delta, *w;
   double alpha_prev, beta, rho_prev, rho;
//   double ztmp;	// [RMS: unused]
   double Theta_prev, Theta, c, sigma_prev, tau_init, tau_prev, tau; 

   /* Parameters used to dynamically update eigenpair */
   double Beta, Delta, Psi, Beta_prev, Delta_prev, Psi_prev;
   double eta;

   double dot_sol, eval_updated, eval_prev, eres2_updated, eres_updated, R;
   double Gamma_prev, Phi_prev;
   double Gamma, Phi;
   double gamma;

   /* The convergence criteria of the inner linear system must satisfy:       */
   /* || current residual || <= relativeTolerance * || initial residual ||    */
   /*                                               + absoluteTol             */

   double relativeTolerance; 
   double absoluteTolerance;
   double LTolerance, ETolerance;

   /* Some constants 							      */
   double tpone = +1.0e+00, tzero = +0.0e+00;

   /* -------------------------------------------*/
   /* Subdivide the workspace into needed arrays */
   /* -------------------------------------------*/

   g      = rwork;
   d      = g + primme->nLocal;
   delta  = d + primme->nLocal;
   w      = delta + primme->nLocal;
   workSpace = w + primme->nLocal;  // This needs at least 2*numOrth+NumEvals)
   
   workSpaceSize = rworkSize - (int)(workSpace - rwork);
   
   /* -----------------------------------------*/
   /* Set up convergence criteria by Tolerance */
   /* -----------------------------------------*/

   if (primme->aNorm <= 0.0L) {
      absoluteTolerance = aNormEstimate*machEps;
      eresTol = eresTol*aNormEstimate;
   }
   else {
      absoluteTolerance = primme->aNorm*machEps;
   }
   tau_prev = tau_init = *rnorm;       /* Assumes zero initial guess */
   LTolerance = eresTol;

   // Andreas: note that eigenresidual tol may not be achievable, because we
   // iterate on P(A-s)P not (A-s). But tau reflects linSys on P(A-s)P.
   if (primme->correctionParams.convTest == primme_adaptive) {
      ETolerance = max(eresTol/1.8L, absoluteTolerance);
      LTolerance = ETolerance;
   }
   else if (primme->correctionParams.convTest == primme_adaptive_ETolerance) {
      LTolerance = max(eresTol/1.8L, absoluteTolerance);
      ETolerance = max(tau_init*0.1L, LTolerance);
   }
   else if (primme->correctionParams.convTest == primme_decreasing_LTolerance) {
      relativeTolerance = pow(primme->correctionParams.relTolBase, 
         (double)-primme->stats.numOuterIterations);
      LTolerance = relativeTolerance * tau_init 
	           + absoluteTolerance + eresTol;
     //printf(" RL %e INI %e abso %e LToler %e aNormEstimate %e \n",
     //relativeTolerance, tau_init, absoluteTolerance,LTolerance,aNormEstimate);
   }
   
   /* --------------------------------------------------------*/
   /* Set up convergence criteria by max number of iterations */
   /* --------------------------------------------------------*/

   /* compute first total number of remaining matvecs */

   maxIterations = primme->maxMatvecs - primme->stats.numMatvecs;

   /* Perform primme.maxInnerIterations, but do not exceed total remaining */
   if (primme->correctionParams.maxInnerIterations > 0) {

      maxIterations = min(primme->correctionParams.maxInnerIterations, 
//.........这里部分代码省略.........

void PithExtractor::process( Billon &billon ) const
{
	billon._pith.clear();

	const int width = billon.n_cols;
	const int height = billon.n_rows;
	const int depth = billon.n_slices;

	uiCoord2D coordPrec, coordCurrent;
	float maxStandList[depth];
	float *maxStandList2 = 0;
	float shift,houghStandThreshold;
	int max,x,y, nbContourPoints;

	max = x = y = 0;

	//extraction de la moelle de la premiere coupe sur la totalité de la coupe
	//nous permet d'avoir un coordonnée a laquelle appliqué la fenetre
	transHough( billon.slice(0), width, height, &x, &y, &max, &nbContourPoints );

	//calcul la coordonnée d'origine de la fenetre
	x = x - (_windowWidth/2);
	x = x < 0 ? 0 : (x > (width - _windowWidth) ? width - _windowWidth : x);
	y = y - (_windowHeight/2);
	y = y < 0 ? 0 : (y > (height - _windowHeight) ? height - _windowHeight : y);
	//extraction de la moelle de la premiere coupe sur la coupe redimensionnée
	coordCurrent = transHough( billon.slice(0), _windowWidth, _windowHeight, &x, &y, &max, &nbContourPoints );

	billon._pith.append(coordCurrent);
	maxStandList[0] = ((float)max)/nbContourPoints;

	//extraction des coupes suivantes
	std::cerr << "extractMoelle : ";
	for(int i=1; i<depth; i++) {
		//~ if(!listeCoupe->contains(i)){
			std::cerr << " " << i;
			coordPrec = billon._pith.last();
			x = x - (_windowWidth/2);
			x = x < 0 ? 0 : (x > (width - _windowWidth) ? width - _windowWidth : x);
			y = y - (_windowHeight/2);
			y = y < 0 ? 0 : (y > (height - _windowHeight) ? height - _windowHeight : y);
			//extraction de la moelle de la coupe i
			coordCurrent = transHough( billon.slice(i), _windowWidth, _windowHeight, &x, &y, &max, &nbContourPoints );
			billon._pith.append(coordCurrent);
			shift = sqrt(pow((double) ((int)(billon._pith.last().x) - (int)(coordPrec.x)),(double) 2.0) + pow( (double)((int)(billon._pith.last().y) - (int)(coordPrec.y)), (double)2.0) );
			//si le resultat obtenu a un decalage trop important avec la coupe precedente alors on recommence l'extraction sur l'ensemble de la coupe
			if(shift > _pithLag && width > _windowWidth && height > _windowHeight){
				std::cerr << "*";
				// std::cerr << "   :> decalage=" << decalage << ">" << _pithLag << "\n";

				billon._pith.pop_back();
				x = y = 0;
				transHough(billon.slice(i), width, height, &x, &y, &max, &nbContourPoints );
				x = x - (_windowWidth/2);
				x = x < 0 ? 0 : (x > (width - _windowWidth) ? width - _windowWidth : x);
				y = y - (_windowHeight/2);
				y = y < 0 ? 0 : (y > (height - _windowHeight) ? height - _windowHeight : y);
				coordCurrent = transHough( billon.slice(i), _windowWidth, _windowHeight, &x, &y, &max, &nbContourPoints );
				billon._pith.append(coordCurrent);
			}
			maxStandList[i] = ((float)max)/nbContourPoints;
			if (!(i % 20)) {
				std::cerr << std::endl;
			}
		//~ }else{
			//~ listMoelle->append(listMoelle->last());
			//~ listMaxStand[i] = 0;
		//~ }
	}
	std::cerr << "\nnbCoupes=" << depth << " listMoelle.size()=" << billon._pith.size() << "\n";

	//calcul du seuil a partir du quel les coupes sont considerées comme erronées
	maxStandList2 = new float[depth];
	memcpy(maxStandList2, maxStandList, depth*sizeof(float));

	qsort(maxStandList2, depth, sizeof(float), &floatCompare);

	int nfc = (_falseCutPercent*depth)/100;
	houghStandThreshold = nfc > 0 ? (maxStandList2[nfc]+maxStandList2[nfc-1])/2 : maxStandList2[nfc];

	delete [] maxStandList2;

	// applique la coorection a la moelle
	correctPith(billon._pith, maxStandList, houghStandThreshold);
}

void Foam::uniformTotalPressureFvPatchScalarField::updateCoeffs
(
    const vectorField& Up
)
{
    if (updated())
    {
        return;
    }

    scalar p0 = pressure_->value(this->db().time().timeOutputValue());

    const fvsPatchField<scalar>& phip =
        patch().lookupPatchField<surfaceScalarField, scalar>(phiName_);

    if (psiName_ == "none" && rhoName_ == "none")
    {
        operator==(p0 - 0.5*(1.0 - pos(phip))*magSqr(Up));
    }
    else if (rhoName_ == "none")
    {
        const fvPatchField<scalar>& psip =
            patch().lookupPatchField<volScalarField, scalar>(psiName_);

        if (gamma_ > 1.0)
        {
            scalar gM1ByG = (gamma_ - 1.0)/gamma_;

            operator==
            (
                p0
               /pow
                (
                    (1.0 + 0.5*psip*gM1ByG*(1.0 - pos(phip))*magSqr(Up)),
                    1.0/gM1ByG
                )
            );
        }
        else
        {
            operator==(p0/(1.0 + 0.5*psip*(1.0 - pos(phip))*magSqr(Up)));
        }
    }
    else if (psiName_ == "none")
    {
        const fvPatchField<scalar>& rho =
            patch().lookupPatchField<volScalarField, scalar>(rhoName_);

        operator==(p0 - 0.5*rho*(1.0 - pos(phip))*magSqr(Up));
    }
    else
    {
        FatalErrorInFunction
            << " rho or psi set inconsitently, rho = " << rhoName_
            << ", psi = " << psiName_ << ".\n"
            << "    Set either rho or psi or neither depending on the "
               "definition of total pressure.\n"
            << "    Set the unused variables to 'none'.\n"
            << "    on patch " << this->patch().name()
            << " of field " << this->dimensionedInternalField().name()
            << " in file " << this->dimensionedInternalField().objectPath()
            << exit(FatalError);
    }

    fixedValueFvPatchScalarField::updateCoeffs();
}

int
main (int argc, char ** argv)
{
  if (argc != 2)
  {
    std::cout << "Usage: " << argv[0] << " [depthMap.yml]" << std::endl;
    return 1;
  }
  IplImage * tmp = (IplImage*)cvLoad(argv[1]);
  cv::Mat depth (tmp);

  cv::Mat boundary (depth.rows, depth.cols, CV_8UC3);

  enum RegionType {Shadow, Veil, Object};

  for (size_t j = 0; j < boundary.rows; ++j)
  {
    for (size_t i = 0; i < boundary.cols; ++i)
    {
      boundary.at<cv::Vec3b>(j, i)[0] = boundary.at<cv::Vec3b>(j, i)[1] = boundary.at<cv::Vec3b>(j, i)[2] = 0;
      if (depth.at<float>(j, i) <= 0)
      {
        boundary.at<cv::Vec3b>(j, i)[Shadow] = 255;
      }
      else
      {
        boundary.at<cv::Vec3b>(j, i)[Object] = 255; 
      }
    }
  }

  cv::Mat mean (boundary.rows, boundary.cols, CV_32FC1);
  cv::Mat var (boundary.rows, boundary.cols, CV_32FC1);

  int window = 2;
  float sumVar = 0;
  int nVar = 0;
  for (size_t j = 0; j < boundary.rows; ++j)
  {
    for (size_t i = 0; i < boundary.cols; ++i)
    {
      if (boundary.at<cv::Vec3b>(j, i)[Shadow] == 255)
      {
        continue;
      }
      float sum = 0;
      float sum2 = 0;
      int n = 0;
      int left = MAX(0, (int)i - window);
      int right = MIN (boundary.cols - 1, (int)i + window);
      int top = MAX (0, (int)j - window);
      int bottom = MIN (boundary.rows - 1, (int)j + window);

      for (int x = left; x < right; ++x)
      {
        for (int y = top; y < bottom; ++y)
        {
          float horz = fabs(depth.at<float>(y, x) - depth.at<float>(y, x + 1));
          float vert = fabs(depth.at<float>(y, x) - depth.at<float>(y + 1, x));
          if (boundary.at<cv::Vec3b>(y, x + 1)[Shadow] == 0)
          {
            sum += horz;
            sum2 += horz*horz;
            ++n;
          }
          if (boundary.at<cv::Vec3b>(y + 1, x)[Shadow] == 0)
          {
            sum += vert;
            sum2 += vert*vert;
            ++n;
          }
        }
      }
     
      if (n)
      {
        mean.at<float>(j, i) = sum/n;
        var.at<float>(j, i) = sum2/n - pow (mean.at<float>(j, i), 2.0);
        sumVar += var.at<float>(j, i);
        ++nVar;
      }
      else
      {
        mean.at<float>(j, i) = 0;
        var.at<float>(j, i) = 0;
      }
    }
  }
  float meanVar = sumVar/nVar;
  for (size_t j = 0; j < boundary.rows; ++j)
  {
    for (size_t i = 0; i < boundary.cols; ++i)
    {
      if (var.at<float>(j, i) > meanVar)
      {
        boundary.at<cv::Vec3b>(j, i)[Veil] = 255;
      }
    }
  }
  cv::imwrite("foo.jpg", boundary);
//.........这里部分代码省略.........

void sdvd(int g[2][K], float es_ovr_n0, long channel_length,
            float *channel_output_vector, int *decoder_output_matrix) {
 
    int i, j, l, ll;                          /* loop variables */
    long t;                                   /* time */
    int memory_contents[K];                   /* input + conv. encoder sr */
    int input[TWOTOTHEM][TWOTOTHEM];          /* maps current/nxt sts to input */
    int output[TWOTOTHEM][2];                 /* gives conv. encoder output */
    int nextstate[TWOTOTHEM][2];              /* for current st, gives nxt given input */
    int accum_err_metric[TWOTOTHEM][2];       /* accumulated error metrics */
    int state_history[TWOTOTHEM][K * 5 + 1];  /* state history table */
    int state_sequence[K * 5 + 1];            /* state sequence list */
 
    int *channel_output_matrix;               /* ptr to input matrix */

    int binary_output[2];                     /* vector to store binary enc output */
    int branch_output[2];                     /* vector to store trial enc output */
 
    int m, n, number_of_states, depth_of_trellis, step, branch_metric,
        sh_ptr, sh_col, x, xx, h, hh, next_state, last_stop; /* misc variables */
 
/* ************************************************************************** */
 
    /* n is 2^1 = 2 for rate 1/2 */
    n = 2;
 
    /* m (memory length) = K - 1 */
    m = K - 1;

    /* number of states = 2^(K - 1) = 2^m for k = 1 */
    number_of_states = (int) pow(2, m);

    /* little degradation in performance achieved by limiting trellis depth
       to K * 5--interesting to experiment with smaller values and measure
       the resulting degradation. */
    depth_of_trellis = K * 5;

    /* initialize data structures */
    for (i = 0; i < number_of_states; i++) {
        for (j = 0; j < number_of_states; j++)
            input[i][j] = 0;

        for (j = 0; j < n; j++) {
            nextstate[i][j] = 0;
            output[i][j] = 0;
        }

        for (j = 0; j <= depth_of_trellis; j++) {
            state_history[i][j] = 0;
        }

        /* initial accum_error_metric[x][0] = zero */
        accum_err_metric[i][0] = 0;
        /* by setting accum_error_metric[x][1] to INT_MAX, we don't need a flag */
        /* so I don't get any more questions about this:  */
        /* INT_MAX is simply the largest possible integer, defined in limits.h */
        accum_err_metric[i][1] = INT_MAX;

    }

    /* generate the state transition matrix, output matrix, and input matrix
       - input matrix shows how FEC encoder bits lead to next state
       - next_state matrix shows next state given current state and input bit
       - output matrix shows FEC encoder output bits given current presumed
       encoder state and encoder input bit--this will be compared to actual
       received symbols to determine metric for corresponding branch of trellis
    */

    for (j = 0; j < number_of_states; j++) {
        for (l = 0; l < n; l++) {
            next_state = nxt_stat(j, l, memory_contents);
            input[j][next_state] = l;

            /* now compute the convolutional encoder output given the current
               state number and the input value */
            branch_output[0] = 0;
            branch_output[1] = 0;

            for (i = 0; i < K; i++) {
                branch_output[0] ^= memory_contents[i] & g[0][i];
                branch_output[1] ^= memory_contents[i] & g[1][i];
            }

            /* next state, given current state and input */
            nextstate[j][l] = next_state;
            /* output in decimal, given current state and input */
            output[j][l] = bin2deci(branch_output, 2);

        } /* end of l for loop */

    } /* end of j for loop */

    #ifdef DEBUG
    printf("\nInput:");
    for (j = 0; j < number_of_states; j++) {
        printf("\n");
        for (l = 0; l < number_of_states; l++)
            printf("%2d ", input[j][l]);
    } /* end j for-loop */

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

static void LIQSS2_recomputeNextTimes(uinteger STATES, LIQSS2_quantizerState p, modelica_real t, modelica_real ft, int dt_index){
    modelica_real diffQ = 0, timeaux = 0;
    uinteger i =0, j =0;
    for(j=0;j<p->nSD[dt_index];j++){
        i= p->SD[dt_index][j];
    // for(j=0;j<STATES;j++){
    //     i=j;
        // if(t>0.00001 && t<0.00001001){
        //     printf("i:  %d\n", i);
        //     printf("t:  %.20f\n", t);
        //     printf("p->nTime[i]:   %.20f \n", p->nTime[i]);
        // }
        if(1){
            // printf("LIQSS2_recomputeNextTimes\n");
            // printf("i:   %d\n", i);
            // printf("t:   %.20f\n", t);
            // printf("p->nTime[i]:   %.20f\n", p->nTime[i]);
            // printf("p->q[i][1]:   %.20f \n", p->q[i][1]);
            // printf("p->x[i][1]:   %.20f \n", p->x[i][1]);
        }
        if(p->ltq[i] == t){
            diffQ = p->q[i][0] - p->qAux[i];
            if(fabs(diffQ) > p->lqu[i] * 1e-6){
                p->a[i] = (p->x[i][1] - p->oldDx[i]) / diffQ;
                if(p->a[i] > 0)
                    p->a[i] = 0;

                // printf("! a[var]: %.20f\n", p->a[i]);
                // printf("! x[cf1]: %.20f\n", p->x[i][1]);
                // printf("! oldDx[var]: %.20f\n", p->oldDx[i]);
                // printf("! diffQ: %.20f\n", diffQ);
            }
        }
        else
            p->flag3[i] = false;

        p->u0[i] = p->x[i][1] - p->q[i][0] * p->a[i];
        p->u1[i] = 2 * p->x[i][2] - p->q[i][1] * p->a[i];
        p->lt[i] = t;
        if(p->flag4[i])
            p->nTime[i] = t;
        else{
            p->diffxq[i][1] = p->q[i][1] - p->x[i][1];
        //     if(1){
            // printf("p->q[i][1] %.20f\n", p->q[i][1]);
            // printf("p->x[i][1] %.20f\n", p->x[i][1]);
        //     printf("p->diffxq[i][1] %.20f\n", p->diffxq[i][1]);
        // }
            p->diffxq[i][2] = -p->x[i][2];
            p->diffxq[i][0] = p->q[i][0] - p->dq[i] + p->lqu[i] - p->x[i][0];
            p->nTime[i] = t + minRootPos(p->diffxq,i,2);
            // printf("nTime[var]: %.20f\n", p->nTime[i]);
            p->diffxq[i][0] = p->q[i][0] - p->dq[i] - p->lqu[i] - p->x[i][0];
            timeaux = t + minRootPos(p->diffxq,i,2);
            // printf("timeaux: %.20f\n", timeaux);
            if(timeaux < p->nTime[i])
                p->nTime[i] = timeaux;
            if(p->a[i] != 0 && fabs(p->x[i][2]) > 1e-10 && !p->flag3[i] && !p->flag2[i]){
                modelica_real coeff[2];
                coeff[0] = p->a[i] * p->a[i] * p->q[i][0] + p->a[i] * p->u0[i] + p->u1[i];
                coeff[1] = p->a[i] * p->a[i] * p->q[i][1] + p->a[i] * p->u1[i];
                // printf("coeff[0] %.20f\n", coeff[0]);
                // printf("coeff[1] %.20f\n", coeff[1]);
                timeaux = t + minRootPosCoeff(&coeff, i, 1);
                // printf("??timeaux: %.20f\n", timeaux);
                if(timeaux < p->nTime[i]){
                    p->flag2[i] = true;
                    p->nTime[i] = timeaux;
                    p->lquOld[i] = p->lqu[i];
                }
            }
            else
                p->flag2[i] = false;

            if(p->nTime[i] > ft)
                p->nTime[i] = ft;

            modelica_real err1 = p->q[i][0] - p->x[i][0] + (p->diffxq[i][1] * (p->nTime[i] - t) / 2) + p->diffxq[i][2] * pow ((p->nTime[i] - t) / 2, 2);
            if (fabs (err1) > 3 * fabs (p->lqu[i])){
                // printf("var: %d\n", i);
                // printf("err1 %.20f\n", err1);

                // printf("t %.20f\n\n", t);
                // printf("lqu[var] %.20f\n", p->lqu[i]);
        	    p->nTime[i] = t + ft * 1.0e-14;
        	}
        }
        // if(t>0.00001 && t<0.00001001){
        //     //printf("Now we are in LIQSS2_recomputeNextTimes\n");
        //     //printf("We are considering the case where i=6\n");
        //     printf("p->nTime[i]:   %.20f \n", p->nTime[i]);
        //     printf("p->diffxq[i][1]:   %.20f \n", p->diffxq[i][1]);
        //     printf("p->diffxq[i][2]:   %.20f \n", p->diffxq[i][2]);
        //     printf("p->diffxq[i][0]:   %.20f \n", p->q[i][0] - p->dq[i] + p->lqu[i] - p->x[i][0]);
        //     printf("p->diffxq[i][0]:   %.20f \n", p->q[i][0] - p->dq[i] - p->lqu[i] - p->x[i][0]);
        // }
        if(1){
            //printf("Now we are in LIQSS2_recomputeNextTimes\n");
            //printf("We are considering the case where i=6\n");
            // printf("p->q[i][0]:   %.20f \n", p->q[i][0]);
//.........这里部分代码省略.........

static double UnquantizeSendRate(guint16 send_rate)
{
	return (send_rate >> 4) * 10.0 / 4096.0 * pow(10.0, (send_rate & 0x000f));
}

double Point::Distance() const //calculating distance from origin as sqrt(x^2+y^2)
{
    return(sqrt(pow(m_x,2)+pow(m_y,2)));
}

double Point::Distance(const Point& Pt) const //calculating distance between two points
{
    return (sqrt(pow(Pt.m_x-m_x,2)+pow(Pt.m_y-m_y,2)));
}

	double p_f(double n) {
		return pow(3.0 * pi * pi * D * n, 1.0 / 3.0);
	}

	double t_E(std::vector<double> n, double f, set_const * C) {
		double res = 0.5 * pow(m_n * m_n * f / C->C_s, 2.0)*C->eta_s(f);

		res += C->U(f);

		int num = n.size();
		/*std::vector<double> nvec;
		for (int i = 0; i < num; i++){
			nvec.push_back(boost::python::extract<double>(n[i]));
		}*/

		double o_sum = 0;
		double r_sum = 0;

		for (int i = 0; i < num; i++){
//			cout << i << ", n[i] = " << n[i] << ", pf = " <<  p_f(n[i]) << endl;
			//res += KineticIntegral(p_f(n[i]), C->M[i] * C->phi_n(C->X_s[i] * f), C);
			res += KineticIntegral(p_f(n[i]), C->M[i] * C->phi_n(C->X_s[i] * (C->M[0]/C->M[i]) * f), C);
			o_sum += C->X_o[i] * n[i];
			r_sum += C->X_r[i] * C->T[i] * n[i];
		}


		res += 0.5 * pow(C->C_o * D * (o_sum) / m_n, 2.0) / (C->eta_o(f));

		res += 0.5 * pow(C->C_r * D * (r_sum) / m_n, 2.0) / (C->eta_r(f));


		//double me = m_e;

		//double pf_e = 0;
		//if (pow(mu_e(n[0] + n[1], n[1], f, C), 2.0) - me*me >= 0){
		//	pf_e = sqrt(pow(mu_e(n[0] + n[1], n[1], f, C), 2.0) - me*me);
		//}

		//if (n[1] != 0.0){
		//	res += KineticIntegral(pf_e, m_e, C);
		//}

		//double mmu = m_mu;
		//double pf_mu = 0;
		//if (pow(mu_e(n[0] + n[1], n[1], f, C), 2.0) - mmu*mmu >= 0){
		//	pf_mu = sqrt(pow(mu_e(n[0] + n[1], n[1], f, C), 2.0) - mmu*mmu);
		//}

		//if (n[1] != 0.0){
		//	res += KineticIntegral(pf_mu, m_mu, C);
		//}

		//gsl_function F;
		//F.function = &func_e;
		//double result = 0.0;
		//double error;
		//double me = m_e;
		//F.params = &me;
		//double pf_e = 0;
		//if (pow(mu_e_2(n[0] + n[1], n[1], f, C), 2.0) - me*me >= 0){
		//	pf_e = sqrt(pow(mu_e_2(n[0] + n[1], n[1], f, C), 2.0) - me*me);
		//}
		////	cout << "MU_E = " <<  mu_e(nn+np, np, f, C) << endl;
		//gsl_integration_qags(&F, 0.0, pf_e, 0, 1e-10, 1000, w_t_E, &result, &error);

		//if (n[1] != 0.0){
		//	//cout << "hey! " << np << endl;
		//	res += result / (pi*pi);
		//}
		////cout << "RESULT   " << result << endl;

		//double mmu = m_mu;
		//F.params = &mmu;
		//double pf_mu = 0;
		//if (pow(mu_e_2(n[0] + n[1], n[1], f, C), 2.0) - mmu*mmu >= 0){
		//	pf_mu = sqrt(pow(mu_e_2(n[0] + n[1], n[1], f, C), 2.0) - mmu*mmu);
		//}
		//gsl_integration_qags(&F, 0.0, pf_mu, 0, 1e-10, 1000, w_t_E, &result, &error);


		//if (n[1] != 0.0){
		//	//cout << "hey! " << np << endl;
		//	res += result / (pi*pi);
		//}

		return res;
	}

ImageObject *
CreateDropShadowMask(const ViewportInfo *vpInfo, int offx, int offy, int r)
{
    int rx2 = r * 2;
    int img_w = vpInfo->file->Width();
    int img_h = vpInfo->file->Height();

    ImageObject *mask = new ImageObject(
        offx + r + img_w, offy + r + img_h, 1);

    // Create the blur image
    float *blur = new float[rx2 * rx2];
    //float root2 = sqrt(2.);
    float blurSum = 0.;
    for (int i=0; i<rx2; i++)
    {
        for (int j=0; j<rx2; j++)
        {
            float u = float(i)/float(rx2-1);
            float v = float(j)/float(rx2-1);

            float gu = exp(-pow((4*(u - .5)), 2));
            float gv = exp(-pow((4*(v - .5)), 2));
            float a = gu * gv;
            blur[i*rx2+j] = a;
            blurSum += a;
        }
    }
    for(int j = 0; j < rx2; ++j)
        for(int i = 0; i < rx2; ++i)
            blur[j*rx2 + i] /= blurSum;

    if(vpInfo->opaqueMode == M_OPAQUE)
    {
        int white = 255;
        mask->SetBlockToColor(offx, offy, offx + img_w, offy + img_h, &white);
    }
    else if(vpInfo->opaqueMode == M_TRANSPARENT)
    {
        int gray = int(vpInfo->opacity * 255.);
        mask->SetBlockToColor(offx, offy, offx + img_w, offy + img_h, &gray);
    }
    else
    {
        unsigned char rr, rg, rb;
        rr = vpInfo->transparentColor[0];
        rg = vpInfo->transparentColor[1];
        rb = vpInfo->transparentColor[2];

        // Everywhere that does not match the bg color, color white.
        for(int y = 0; y < img_h; ++y)
        {
            for(int x = 0; x < img_w; ++x)
            {
                unsigned char *src = vpInfo->file->Pixel(x, y);
                unsigned char *dest = mask->Pixel(x+offx, y+offy);

                if(src[0] != rr || src[1] != rg || src[2] != rb)
                    *dest = 255;
            }
        }
    }

    // Blur the mask by convolving it with the blur kernel.
    for(int y = offy; y < offy + img_h; ++y)
    {
        for(int x = offx; x < offx + img_w; ++x)
        {
            int ksrcx = x - r;
            int ksrcy = y - r;
            float channelSum = 0.;
            for(int ky = 0; ky < rx2; ++ky)
            {
                for(int kx = 0; kx < rx2; ++kx)
                {
                    unsigned char *mask_ptr = mask->Pixel(ksrcx + kx, ksrcy + ky);
                    channelSum += float(*mask_ptr) * blur[ky*rx2+kx];
                }
            }
            unsigned char *dest = mask->Pixel(x, y);
            *dest = (unsigned char)(int)(channelSum);
        }
    }

    if(vpInfo->opaqueMode == M_TRANSPARENT)
    {
        int black = 0;
        mask->SetBlockToColor(0, 0, img_w, img_h, &black);
    }
    delete [] blur;

    return mask;
}

//.........这里部分代码省略.........
			//float prevphase = fmod(unit->m_prevphase[index]+phaseadvance,TWOPI);
			//float a,b,tmp;
			//a=phasenow; b=prevphase;
			//b=phasenow; a=prevphase;

			//if(b<a) {b= b+TWOPI;}

			//float phasechange = sc_min(b-a,a+TWOPI-b); //more complicated, need mod 2pi and to know lower and upper
			//phasesum+= phasechange;
			//unit->m_prevphase[index]= phasenow;

			//((p->bin[index-1].mag) * (p->bin[index-1].mag))

			//printf("comparison %f %f \n",fftbuf[g_bins2[index]], power);
			//sum+= (unit->m_weights[index])* power;

			sum+= (weights[index])* (fftbuf[bins[index]]);
		}


		//transient test here too?
		//if(phasesum>(5*PI)){sum=0.0;}

		//if((i>5) && (i<15))
		//printf("test phasesum %f \n", phasesum);
		//unit->m_leaknote[i] = (0.8*unit->m_leaknote[i]) + sum;

		chroma[chromaindex]+= sum; //unit->m_leaknote[i]; //sum;
	}

	float* key = unit->m_key;

	//major
	for (int i=0;i<12;++i) {

		sum=0.0;
		for (int j=0;j<7;++j) {
			indexbase=g_major[j];

			index=(i+indexbase)%12;
			//sum+=(chroma[index]*g_kkmajor[indexbase]);

			sum+=(chroma[index]*g_diatonicmajor[indexbase]);

		}

		key[i]=sum; //10*log10(sum+1);
	}

	//minor
	for (int i=0;i<12;++i) {

		sum=0.0;
		for (int j=0;j<7;++j) {
			indexbase=g_minor[j];

			index=(i+indexbase)%12;
			//sum+=(chroma[index]*g_kkminor[indexbase]);

			sum+=(chroma[index]*g_diatonicminor[indexbase]);

		}

		key[12+i]=sum;
	}

	float keyleak= ZIN0(1); //fade parameter to 0.01 for histogram in seconds, convert to FFT frames

	//keyleak in seconds, convert to drop time in FFT hop frames (FRAMEPERIOD)
	keyleak= sc_max(0.001f,keyleak/unit->m_frameperiod); //FRAMEPERIOD;

	//now number of frames, actual leak param is decay exponent to reach 0.01 in x seconds, ie 0.01 = leakparam ** (x/ffthopsize)
	//0.01 is -40dB
	keyleak= pow(0.01f,(1.f/keyleak));

	float * histogram= unit->m_histogram;

	int bestkey=0;
	float bestscore=0.0;

	for (int i=0;i<24;++i) {
		histogram[i]= (keyleak*histogram[i])+key[i];

		if(histogram[i]>bestscore) {
			bestscore=histogram[i];
			bestkey=i;
		}

	//printf("%f ",histogram[i]);
	}

	//should find secondbest and only swap if win by a margin

	//printf(" best %d \n\n",bestkey);
	//what is winning currently? find max in histogram
	unit->m_currentKey=bestkey;

	//about 5 times per second
	//if((unit->m_triggerid) && ((unit->m_frame%2==0))) SendTrigger(&unit->mParent->mNode, unit->m_triggerid, bestkey);
}

//.........这里部分代码省略.........
                    for( int m=1; m <= mpts; m++ )
                    {

                        double a   = spts.get(m) - 2.0 * large_eta;
                        double b   = spts.get(m);

                        int ishift = get_shift_index( a );

                        double integral = 0.0;
                        if( ishift == 0 )
                        {
                            // no shift takes place, can simply integrate from a to b
                            for( int kq = 1; kq <= kmax; kq++ )
                            {
                                integral += qold.get(i, me, kq) * 
                                    dgc->integratePhi( a, b, kq );
                            }

                        }
                        else if( ishift < 0 )
                        {

                            // integral at the shifted value
                            for( int kq = 1; kq <= kmax; kq++ )
                            {
                                integral += qold.get(i+ishift, me, kq) * 
                                    dgc->integratePhi( shift_xi(a), 1.0, kq );
                            }

                            // integrals between i and ishift
                            for( int n=-1; n > ishift; n-- )
                            { integral += qold.get(i + n, meqn, 1) * dgc->integratePhi(-1.0, 1.0, 1); }

                            // integral inside current cell
                            for( int kq = 1; kq <= kmax; kq++ )
                            {
                                integral += qold.get(i, me, kq) * 
                                    dgc->integratePhi( -1.0, spts.get(m), kq );
                            }

                        }
                        else
                        {
                            printf("   StepAdvecFluxForm.cpp problem can't handle negative velocities yet\n");
                            exit(1);
                        }

                        qnp1 += 0.25 * dx * wgt.get(m) * phi_x.get(m,k) * integral;

                    }
                }

                // left flux
                double a = -1.0 - 2.0 * large_eta;
                int ishift = get_shift_index( a );

                double Fm = 0.0;
                if( ishift < 0 )
                {
                    for( int kq = 1; kq <= kmax; kq++ )
                    {
                        Fm += qold.get(i+ishift, me, kq) * 
                            dgc->integratePhi( shift_xi(a), 1.0, kq ) / 2.0;
                    }
                    for( int n=-1; n > ishift; n-- )
                    {
                        Fm += qold.get(i+n, me, 1);
                    }
                }
                Fm *= sqrt(2*k-1) * pow(-1, k+1 );

                // right flux
                double Fp = 0.0;
                a += 2.0;
                ishift = get_shift_index( a );
                if( ishift < 1 )
                {
                    for( int kq = 1; kq <= kmax; kq++ )
                    {
                        Fp += qold.get(i+ishift, me, kq) * 
                            dgc->integratePhi( shift_xi(a), 1.0, kq ) / 2.0;
                    }
                }
                for( int n=0; n > ishift; n-- )
                { Fp += qold.get(i+n, me, 1); }
                Fp *= sqrt(2*k-1);

                qnew.set(i, me, k, qnp1 - ( Fp - Fm ) );

            }

            delete dgc;

        }//end of loop over each cell
    }//end of loop over each eqn  

    void SetBndValues(const dTensor2& node, dTensorBC3& aux, dTensorBC3& q);
    SetBndValues(node, auxvals, qnew);

}

body_part Creature::select_body_part(Creature *source, int hit_roll)
{
    // Get size difference (-1,0,1);
    int szdif = source->get_size() - get_size();
    if(szdif < -1) {
        szdif = -1;
    } else if (szdif > 1) {
        szdif = 1;
    }

    if(g->debugmon) {
        g->add_msg("source size = %d", source->get_size());
        g->add_msg("target size = %d", get_size());
        g->add_msg("difference = %d", szdif);
    }

    std::map<body_part, double> hit_weights = default_hit_weights[szdif];
    std::map<body_part, double>::iterator iter;

    // If the target is on the ground, even small/tiny creatures may target eyes/head. Also increases chances of larger creatures.
    // Any hit modifiers to locations should go here. (Tags, attack style, etc)
    if(is_on_ground()) {
        hit_weights[bp_eyes] += 10;
        hit_weights[bp_head] += 20;
    }

    //Adjust based on hit roll: Eyes, Head & Torso get higher, while Arms and Legs get lower.
    //This should eventually be replaced with targeted attacks and this being miss chances.
    hit_weights[bp_eyes] = floor(hit_weights[bp_eyes] * pow(hit_roll, 1.15) * 10);
    hit_weights[bp_head] = floor(hit_weights[bp_head] * pow(hit_roll, 1.15) * 10);
    hit_weights[bp_torso] = floor(hit_weights[bp_torso] * pow(hit_roll, 1) * 10);
    hit_weights[bp_arms] = floor(hit_weights[bp_arms] * pow(hit_roll, 0.95) * 10);
    hit_weights[bp_legs] = floor(hit_weights[bp_legs] * pow(hit_roll, 0.975) * 10);


    // Debug for seeing weights.
    if(g->debugmon) {
        g->add_msg("eyes = %f", hit_weights.at(bp_eyes));
        g->add_msg("head = %f", hit_weights.at(bp_head));
        g->add_msg("torso = %f", hit_weights.at(bp_torso));
        g->add_msg("arms = %f", hit_weights.at(bp_arms));
        g->add_msg("legs = %f", hit_weights.at(bp_legs));
    }

    double totalWeight = 0;
    std::set<weight_pair, weight_compare> adjusted_weights;
    for(iter = hit_weights.begin(); iter != hit_weights.end(); ++iter) {
        totalWeight += iter->second;
        adjusted_weights.insert(*iter);
    }

    double roll = rng_float(1, totalWeight);
    body_part selected_part = bp_torso;

    std::set<weight_pair, weight_compare>::iterator adj_iter;
    for(adj_iter = adjusted_weights.begin(); adj_iter != adjusted_weights.end(); ++adj_iter) {
        roll -= adj_iter->second;
        if(roll <= 0) {
            selected_part = adj_iter->first;
            break;
        }
    }

    return selected_part;
}

//---------------------------------------------------------
bool CWatersheds_ext::Get_Basin(CSG_Grid *pBasins, CSG_Shapes *pPolygons, int xMouth, int yMouth, int Main_ID)
{
	int						x, y, Basin_ID	= 1 + pPolygons->Get_Count();
	CSG_Shape				*pPolygon;
	CSG_Grid_Stack			Stack;
	CSG_Simple_Statistics	s_Height, s_Distance;

	//-----------------------------------------------------
	Stack.Push(x = xMouth, y = yMouth);

	pBasins		->Set_Value(x, y, Basin_ID);
	m_Distance	 .Set_Value(x, y, 0.0);

	s_Height	+= m_pDEM->asDouble(x, y);
	s_Distance	+= 0.0;

	//-----------------------------------------------------
	while( Stack.Get_Size() > 0 && Process_Get_Okay() )
	{
		Stack.Pop(x, y);

		double	d	= m_Distance.asDouble(x, y);

		//-------------------------------------------------
		for(int i=0; i<8; i++)
		{
			int	ix	= Get_xFrom(i, x);
			int	iy	= Get_yFrom(i, y);

			if( is_InGrid(ix, iy) && pBasins->is_NoData(ix, iy) && i == m_Direction.asInt(ix, iy) )
			{
				Stack.Push(ix, iy);

				pBasins		->Set_Value(ix, iy, Basin_ID);
				m_Distance	 .Set_Value(ix, iy, d + Get_Length(i));

				s_Height	+= m_pDEM->asDouble(ix, iy);
				s_Distance	+= d + Get_Length(i);
			}
		}
	}

	//-----------------------------------------------------
	if( s_Height.Get_Count() > 1 && (pPolygon = Get_Basin(pBasins, pPolygons)) != NULL )
	{
		double		d, Area, Perimeter, Side_A, Side_B;
		CSG_String	Gravelius;

	//	Area		= s_Height.Get_Count() * Get_System()->Get_Cellarea();
		Area		= ((CSG_Shape_Polygon*)pPolygon)->Get_Area();
		Perimeter	= ((CSG_Shape_Polygon*)pPolygon)->Get_Perimeter();

		d			= 0.28 * Perimeter / sqrt(Area);
		Gravelius	= d > 1.75 ? _TL("rectangular")
					: d > 1.5  ? _TL("ovalooblonga-rectangularoblonga")
					: d > 1.25 ? _TL("ovaloredonda-ovalooblonga")
					:            _TL("redonda-ovaloredonda");

		d			= pow(Perimeter, 2.0) - 8.0 * Area;
		Side_A		= d > 0.0 ? (Perimeter + sqrt(d))      / 4.0 : -1.0;
		Side_B		= d > 0.0 ? (Perimeter - 2.0 * Side_A) / 2.0 : -1.0;

		pPolygon->Set_Value(FIELD_ID			, Basin_ID);
		pPolygon->Set_Value(FIELD_ID_MAIN		, Main_ID);

		pPolygon->Set_Value(FIELD_MOUTH_X		, Get_System()->Get_xGrid_to_World(xMouth));
		pPolygon->Set_Value(FIELD_MOUTH_Y		, Get_System()->Get_yGrid_to_World(yMouth));

		pPolygon->Set_Value(FIELD_PERIMETER		, Perimeter);
		pPolygon->Set_Value(FIELD_AREA			, Area);

		pPolygon->Set_Value(FIELD_CENTROID_X	, ((CSG_Shape_Polygon*)pPolygon)->Get_Centroid().x);
		pPolygon->Set_Value(FIELD_CENTROID_Y	, ((CSG_Shape_Polygon*)pPolygon)->Get_Centroid().y);

		pPolygon->Set_Value(FIELD_Z_MEAN		, s_Height  .Get_Mean());
		pPolygon->Set_Value(FIELD_Z_RANGE		, s_Height  .Get_Range());

		pPolygon->Set_Value(FIELD_DIST_MEAN		, s_Distance.Get_Mean());
		pPolygon->Set_Value(FIELD_DIST_MAX		, s_Distance.Get_Maximum());

		pPolygon->Set_Value(FIELD_CONCTIME		, s_Height.Get_Range() <= 0.0 ? -1.0 :
			pow(0.87 * pow(s_Distance.Get_Maximum() / 1000.0, 3.0) / s_Height.Get_Range(),  0.385)
		);

		pPolygon->Set_Value(FIELD_BASIN_TYPE	, Gravelius);

		pPolygon->Set_Value(FIELD_EQVRECT_A		, Side_A);
		pPolygon->Set_Value(FIELD_EQVRECT_B		, Side_B);

		pPolygon->Set_Value(FIELD_OROG_IDX		, SG_Get_Square(s_Height.Get_Mean()) / (0.0001 * Area));	// Orographic index, defined as the mean catchment altitude times the ratio of the mean catchment altitude to the orthogonal projection of drainage area (Alcázar, Palau (2010): Establishing environmental flow regimes in a Mediterranean watershed based on a regional classification. Journal of Hydrology, V. 388
		pPolygon->Set_Value(FIELD_MASS_IDX		, Perimeter / (0.0001 * Area));								// Perimeter / (0.0001 * Area) ??!!

		pPolygon->Set_Value(FIELD_BASINS_UP		, 0.0);	// Upslope Basins
		pPolygon->Set_Value(FIELD_BASINS_DOWN	, 0.0);	// Downslope Basins

		return( true );
	}

	return( false );
//.........这里部分代码省略.........