C++ VecDoub类代码示例

VecDoub Actions_AxisymmetricFudge_InterpTables::actions(const VecDoub &XX, void* with_f){
    // THIS IS A HORRIBLE FUDGE TO STOP J_R=0
    VecDoub X = XX;
    if(fabs(X[3])<0.001)X[3]=0.001;

    bool wf;
    if(with_f)
        wf=(bool *)with_f;

    double E = Pot->H(X),Lz = Pot->Lz(X),JR,Jz;
    double Delta = UV->findDelta_interp(E,fabs(Lz));
    VecDoub Iuv = IuIv(X,Delta,E,Lz*Lz);
    if(int_acts(Lz,E,Iuv[0],Iuv[1],&JR,&Jz)==0 or no_table){
        double alpha = -1.-Delta*Delta;
        VecDoub JJ = UV->actions(X,&alpha);
        JR=JJ[0]; Jz=JJ[2];
    }
    VecDoub JJ = {JR,Lz,Jz};
    if(wf){
        if(Lz<Lcgrid[0]) JJ.push_back(Pot->R_L(Lz,Rgrid[0]));
        else if(Lz>Lcgrid[NR-1]) JJ.push_back(Pot->R_L(Lz,Rgrid[NR-1]));
        else{
            int bot,top;
            topbottom(Lcgrid, Lz, &bot, &top,"actions");
            JJ.push_back(Rgrid[bot]+(Lz-Lcgrid[bot])/(Lcgrid[top]-Lcgrid[bot])*(Rgrid[top]-Rgrid[bot]));
        }
        VecDoub freqs = PotentialFreq(JJ[3]);
        for(auto i:freqs)JJ.push_back(i);
    }
    return JJ;
}

// ======================================================================================
// Galactic <==> Cartesian
VecDoub GalacticToCartesian(const VecDoub &Galactic,
								      const VecDoub& SolarPosition){
	// l,b,s->X,Y,Z
	double 	cl = cos(Galactic[0]), sl = sin(Galactic[0]),
			cb = cos(Galactic[1]), sb = sin(Galactic[1]);

	double x = Galactic[2]*cb*cl;
	double z = Galactic[2]*sb;
	// Need to rotate to account for the height of the Sun above the plane
	double h = sqrt(SolarPosition[0]*SolarPosition[0]
	                +SolarPosition[1]*SolarPosition[1]);
	double ct = SolarPosition[0]/h, st = SolarPosition[1]/h;

	VecDoub Cartesian {	SolarPosition[0]-ct*x-st*z,
						-Galactic[2]*cb*sl,
						-st*x+ct*z+SolarPosition[1]};
	if(Galactic.size()==3)return Cartesian;
	// vlos,mu_lcos(b),mu_b -> vx,vy,vz
	// in units km/s, mas/yr -> km/s
	else{
		double vl = PM_Const*Galactic[2]*Galactic[4];
		double vb = PM_Const*Galactic[2]*Galactic[5];
		double tmp = cb*Galactic[3]-sb*vb;

		double vx = cl*tmp-sl*vl+SolarPosition[2];
		double vy = sl*tmp+cl*vl+SolarPosition[3];
		double vz = sb*Galactic[3]+cb*vb+SolarPosition[4];
		VecDoub CartVel{-(vx*ct+vz*st),-vy,-vx*st+vz*ct};
	  	for (	VecDoub::iterator it = CartVel.begin();
				it != CartVel.end(); ++it) Cartesian.push_back(*it);
			return Cartesian;
	}
}

VecDoub CartesianToSphericalPolar(const VecDoub& Cartesian){
	double r = sqrt(Cartesian[0]*Cartesian[0]+Cartesian[1]*Cartesian[1]+Cartesian[2]*Cartesian[2]);
    VecDoub SPolar = {r,atan2(Cartesian[1],Cartesian[0]),acos(Cartesian[2]/r)};
    if(Cartesian.size()==3)	return SPolar;
    SPolar.push_back((Cartesian[3]*cos(SPolar[1])+Cartesian[4]*sin(SPolar[1]))*sin(SPolar[2])+cos(SPolar[2])*Cartesian[5]);
    SPolar.push_back(-Cartesian[3]*sin(SPolar[1])+Cartesian[4]*cos(SPolar[1]));
    SPolar.push_back((Cartesian[3]*cos(SPolar[1])+Cartesian[4]*sin(SPolar[1]))*cos(SPolar[2])-sin(SPolar[2])*Cartesian[5]);
    return SPolar;
}

/*! Draw connected line segments in page coordinates (pts), with options to close
 * /and/or fill the curve, or to set the curve as a clip area.
 */
void PSpage::polyline(VecDoub &x, VecDoub &y, Int close, Int fill, Int clip)
{
	Int i,n = min(x.size(),y.size());
	fprintf(PSpage::PLT,"np %g %g mt\n",x[0],y[0]);
	for (i=1;i<n;i++) fprintf(PSpage::PLT,"%g %g lt\n",x[i],y[i]);
	if (close || fill || clip) fprintf(PSpage::PLT,"cp ");
	if (fill) fprintf(PSpage::PLT,"fi\n");
	else if (clip) fprintf(PSpage::PLT,"clip\n");
	else fprintf(PSpage::PLT,"st\n");
}

VecDoub Actions_Spherical::find_limits(double r, double E, double L){
    VecDoub limits;
    Actions_Spherical_limits_struct Act(Pot,E,L);
    double r_in=r, r_out=r;
    root_find RF(SMALL,100);
    if(p_r(0.,&Act)>0) r_in=0.;
    else while(p_r(r_in,&Act)>=0.0) r_in*=0.9;
    while(p_r(r_out,&Act)>=0.0) r_out*=1.1;
    limits.push_back(RF.findroot(&p_r,r_in,r,&Act));
    limits.push_back(RF.findroot(&p_r,r,r_out,&Act));
    return limits;
}

VecDoub StackelTriaxial::tau2ints(const VecDoub& tau){
	VecDoub pp = CS->tau2p(tau);
	double X = 0.5*pp[0]-(tau[0]+CS->alpha())*(tau[0]+CS->gamma())*G(tau[0])/(tau[0]-tau[1])/(tau[0]-tau[2]);
	double Y = 0.5*pp[1]-(tau[1]+CS->alpha())*(tau[1]+CS->gamma())*G(tau[1])/(tau[1]-tau[0])/(tau[1]-tau[2]);
	double Z = 0.5*pp[2]-(tau[2]+CS->alpha())*(tau[2]+CS->gamma())*G(tau[2])/(tau[2]-tau[1])/(tau[2]-tau[0]);
	VecDoub Ints = {X+Y+Z};
	double J =(tau[1]+tau[2])*X+(tau[2]+tau[0])*Y+(tau[0]+tau[1])*Z;
	double K = tau[1]*tau[2]*X+tau[2]*tau[0]*Y+tau[0]*tau[1]*Z;
	Ints.push_back((CS->alpha()*(CS->alpha()*Ints[0]+J)+K)/(CS->alpha()-CS->gamma()));
	Ints.push_back((CS->gamma()*(CS->gamma()*Ints[0]+J)+K)/(CS->gamma()-CS->alpha()));
	return Ints;
}

void Eigsym::eig(const MatDoub &A, MatDoub &V, VecDoub &lambda) {
	unsigned int n = A.ncols();	/* size of the matrix */
	double a[n*n]; 	/* store initial matrix */
	double w[n];		/* store eigenvalues */
	int matz = 1;		/* return both eigenvalues as well as eigenvectors */
	double x[n*n]; 	/* store eigenvectors */

	for(unsigned int i=0; i<n; i++) {
		for(unsigned int j=0; j<n; j++) {
			a[i*n+j] = A[i][j];
		}
	}

	unsigned int ierr = 0;

	ierr = rs ( n, a, w, matz, x );

	V.assign(n,n,0.0);
	lambda.resize(n);
	for(unsigned int i=0; i<n; i++) {
		lambda[i] = w[i];
		for(unsigned int j=0; j<n; j++) {
			V[j][i] = x[i*n+j];
		}
	}
}

VecDoub SphericalPolarToCartesian(const VecDoub& Spherical){
    VecDoub SPolar = {Spherical[0]*sin(Spherical[2])*cos(Spherical[1]),
    					Spherical[0]*sin(Spherical[2])*sin(Spherical[1]),
    					Spherical[0]*cos(Spherical[2])};
    if(Spherical.size()==3)	return SPolar;
    return SPolar;
}

std::vector<VecDoub> EquatorialToGalacticwithErrors(const VecDoub &Equatorial,
	const VecDoub &Errors){
	//alpha, dec, s => l,b,s
	double alpha = Equatorial[0], delta = Equatorial[1];
	double cd = cos(delta), sd = sin(delta);
	double dalpha = alpha-RA_GP;
	double b=asin(sdGP*sd+cdGP*cd*cos(dalpha));
	double l=lCP-atan2(cd*sin(alpha-RA_GP),cdGP*sd-sdGP*cd*cos(dalpha));
	if(l<0.)l+=2.*PI;
	if(Equatorial.size()==3){	std::vector<VecDoub> Galactic {{l,b,Equatorial[2]},Errors};
								return Galactic;}
	else{
		//vlos, ma_cos(d), md => vlos, ml_cos(b), mb
		double cb = cos(b), sb = sin(b);
		double A11=(sdGP*cd-cdGP*sd*cos(dalpha))/cb;
		double A12=-cdGP*sin(dalpha)/cb;
		double A21,A22;
		double dl = lCP-l;
		if(fabs(cos(dl))>fabs(sin(dl))){
			A21=(sd*sin(dalpha)-sb*sin(dl)*A11)/cos(dl);
			A22=-(cos(dalpha)+sb*sin(dl)*A12)/cos(dl);
		}else{
			A21=(cdGP*cd+sdGP*sd*cos(dalpha)+sb*cos(dl)*A11)/sin(dl);
			A22=(sdGP*sin(dalpha)+sb*cos(dl)*A12)/sin(dl);
		}

		std::vector<VecDoub> Galactic {
		{l,b,Equatorial[2],Equatorial[3],
		A21*Equatorial[5]+A22*Equatorial[4],A11*Equatorial[5]+A12*Equatorial[4]}
		,{Errors[0],Errors[1],Errors[2],Errors[3],
		sqrt(A21*A21*Errors[5]*Errors[5]+A22*A22*Errors[4]*Errors[4]),
		sqrt(A11*A11*Errors[5]*Errors[5]+A12*A12*Errors[4]*Errors[4])}};
		return Galactic;
		}
}

/*
 Utility method used by the Spectral method fine-tune an initial
 given community split.
 */
void modifySplit( double tol, int countmax ) {

    double qmax  = 0;
    double qold  = 0;
    int count    = 0;
    int Ng       = si.size();

    visited.resize(Ng);

    VecDoub Gsi(Ng);
    MatDoub GSI(Ng,2);

    for(int i=0; i<Ng; i++) {
        Gsi[i]    = si[i];
        GSI[i][0] = SI[i][0];
        GSI[i][1] = SI[i][1];
    }

    maxModularity( qmax );

    while( count < countmax ) {

        if( qmax > qold ) {
            for(int i=0; i<Ng; i++) {
                Gsi[i]    = si[i];
                GSI[i][0] = SI[i][0];
                GSI[i][1] = SI[i][1];
            }
        }

        qold = qmax;
        qmax = 0.0;

        maxModularity(qmax);

        count++;

    }

    for(int i=0; i<Ng; i++) {
        si[i]    = Gsi[i];
        SI[i][0] = GSI[i][0];
        SI[i][1] = GSI[i][1];
    }

}

/*
 Utility method used by the Spectral method to find
 which node, when moved gives the maximum change in the
 Modularity value.
 */
void maxModularity(double &qmax) {

    int N = si.size();

    VecDoub qstored(N);
    for(int i=0; i<N; i++)
        qstored[i]     = 0;


    double Q = 0.0;

    for(int k=0; k<N; k++) {

        if( visited[k] < 1 ) {

            Q          = 0.0;

            deltaModularityMax( k, Q );

            qstored[k] = Q;

        }
    }

    qmax        =  0;//qstored(0);
    int ind_max =  -1;//0;
    for(int i=0; i<N; i++) {

        if( qstored[i] > qmax ) {
            qmax    = qstored[i];
            ind_max = i;
        }

    }


    for(int i=0; i<N; i++) {

        if( i != ind_max )
            ;
        else {
            visited[i] = 1;
            if( si[i] == 1 ) {
                si[i] = -1;
                SI[i][0] = 0;
                SI[i][1] = 1;
            } else {
                si[i] = 1;
                SI[i][0] = 1;
                SI[i][1] = 0;
            }
        }
    }


}

// ============================================================================
// Dehnen Potential
// ============================================================================
double Dehnen::Phi(const VecDoub &x){
	/* potential at Cartesian x */
	assert(x.size()==3);
	double r = norm<double>(x);
	double chi = pow(r/rs,1./alpha);
	chi=chi/(1+chi);
	return -conv::FPG*rhoS*rs*rs*alpha*
	(rs/r*incomplete_beta(alpha*(3-gamma),alpha*(beta-3),chi)
	+incomplete_beta(alpha*(beta-2),alpha*(2-gamma),1-chi));
}

VecDoub GalacticToEquatorial(const VecDoub &Galactic){
	//l,b,s => alpha, dec, s
	double l = Galactic[0], b = Galactic[1];
	double cb = cos(b),sb = sin(b);
	double dl = lCP-l;
	double delta=asin(cdGP*cb*cos(-dl)+sb*sdGP);
	double alpha=RA_GP+atan2(cb*sin(dl),sb*cdGP-cb*sdGP*cos(-dl));
	if(alpha>2.*PI)alpha-=2.*PI;
	VecDoub Equatorial {alpha,delta,Galactic[2]};
	if(Galactic.size()==3)return Equatorial;
	else{
		double dalpha = alpha-RA_GP;
		//vlos, ml_cos(b), mb => vlos, ma_cos(d), md
		double cd = cos(delta), sd = sin(delta);
		double A11=(sdGP*cd-cdGP*sd*cos(dalpha))/cb;
		double A12=-cdGP*sin(dalpha)/cb;
		double A21,A22;
		if(fabs(cos(dl))>fabs(sin(dl))){
			A21=(sd*sin(dalpha)-sb*sin(dl)*A11)/cos(dl);
			A22=-(cos(dalpha)+sb*sin(dl)*A12)/cos(dl);
		}else{
			A21=(cdGP*cd+sdGP*sd*cos(dalpha)+sb*cos(dl)*A11)/sin(dl);
			A22=(sdGP*sin(dalpha)+sb*cos(dl)*A12)/sin(dl);
		}
		double Prod = A11*A22-A12*A21;
		VecDoub EqVel {Galactic[3],(A11*Galactic[4]-A21*Galactic[5])/Prod,
						(A22*Galactic[5]-A12*Galactic[4])/Prod};
		for (	VecDoub::iterator it = EqVel.begin();
				it != EqVel.end(); ++it) Equatorial.push_back(*it);
		return Equatorial;
		}
}

VecDoub EquatorialToGalactic(const VecDoub &Equatorial){
	//alpha, dec, s => l,b,s
	double alpha = Equatorial[0], delta = Equatorial[1];
	double cd = cos(delta), sd = sin(delta);
	double dalpha = alpha-RA_GP;
	double b=asin(sdGP*sd+cdGP*cd*cos(dalpha));
	double l=lCP-atan2(cd*sin(alpha-RA_GP),cdGP*sd-sdGP*cd*cos(dalpha));
	if(l<0.)l+=2.*PI;
	VecDoub Galactic {l,b,Equatorial[2]};
	if(Equatorial.size()==3)return Galactic;
	else{
		//vlos, ma_cos(d), md => vlos, ml_cos(b), mb
		double cb = cos(b), sb = sin(b);
		double dl = lCP-l;
		double A11=(sdGP*cd-cdGP*sd*cos(dalpha))/cb;
		double A12=-cdGP*sin(dalpha)/cb;
		double A21,A22;
		if(fabs(cos(dl))>fabs(sin(dl))){
			A21= (sd*sin(dalpha)-sb*sin(dl)*A11)/cos(dl);
			A22=-(   cos(dalpha)+sb*sin(dl)*A12)/cos(dl);
		}else{
			A21=(cdGP*cd+sdGP*sd*cos(dalpha)+sb*cos(dl)*A11)/sin(dl);
			A22=(sdGP*sin(dalpha)+sb*cos(dl)*A12)/sin(dl);
		}

		VecDoub GalVel {Equatorial[3],A21*Equatorial[5]+A22*Equatorial[4],
						A11*Equatorial[5]+A12*Equatorial[4]};
		for (	VecDoub::iterator it = GalVel.begin();
				it != GalVel.end(); ++it) Galactic.push_back(*it);
		return Galactic;
		}
}

VecDoub CartesianToGalactic(const VecDoub &Cartesian,
									const VecDoub& SolarPosition){
	// X,Y,Z->l,b,s
	double tmp1 = SolarPosition[0]-Cartesian[0];
	double tmp2 = -Cartesian[1];
	double tmp3 = Cartesian[2]-SolarPosition[1];
	// Need to rotate to account for the height of the Sun above the plane
	double h = sqrt(SolarPosition[0]*SolarPosition[0]
	                +SolarPosition[1]*SolarPosition[1]);
	double ct = SolarPosition[0]/h, st = SolarPosition[1]/h;

	double x = tmp1*ct-tmp3*st, z = tmp1*st+tmp3*ct;

	double Distance = norm<double>({x,tmp2,z});

	VecDoub Galactic {	atan2(tmp2,x),
						asin(z/Distance),
						Distance};
	if(Cartesian.size()==3)return Galactic;
	// vx,vy,vz -> vlos,mu_lcos(b),mu_b
	// in units km/s -> km/s mas/yr
	else{ 	double vx=-Cartesian[3]*ct-Cartesian[5]*st-SolarPosition[2];
			double vy = -Cartesian[4]-SolarPosition[3];
			double vz = Cartesian[5]*ct+Cartesian[3]*st-SolarPosition[4];
			double 	cl = cos(Galactic[0]), sl = sin(Galactic[0]),
			cb = cos(Galactic[1]), sb = sin(Galactic[1]);
			VecDoub GalVel {vx*cl*cb+vy*sl*cb+vz*sb,(-vx*sl+vy*cl)/(PM_Const*Distance),
				        	(-vx*cl*sb-vy*sl*sb+vz*cb)/(PM_Const*Distance)};
			for (	VecDoub::iterator it = GalVel.begin();
					it != GalVel.end(); ++it) Galactic.push_back(*it);
			return Galactic;
		}
}