C++ energy函数代码示例

int main(int,char** argv)
{
 
  auto solar_system = construct_tuple(sun,jupiter,saturn,uranus,neptune);
  offset(solar_system);

  printf ("%.9f\n", energy(solar_system));
  
  int n = atoi(argv[1]);
  
  for (int i = 1; i <= n; i++)
  {
    advance(solar_system);
  }

  printf ("%.9f\n", energy(solar_system));

  return 0;
}

void heartbeat(struct reb_simulation* r){
	if (r->t > next_output){
		next_output *= 1.02;
		reb_integrator_synchronize(r);
		FILE* f = fopen("energy.txt","a");
		double e = energy(r);
		fprintf(f, "%e %.16e\n", r->t, fabs((e-e_init)/e_init));
		fclose(f);
	}
}

 EnergyGradAndPrecond energy_grad_and_precond(const Vector &x) const {
   EnergyGradAndPrecond egpg;
   egpg.energy = energy(x);
   egpg.grad = grad(x);
   egpg.precond = grad(x);
   const int sz = x.get_size();
   for (int i=0; i<sz; i++) {
     egpg.precond[i] /= spring(i);
   }
   return egpg;
 }

void problem_output(){
	if (output_check(10000.)){
		output_timing();
		integrator_synchronize();
		FILE* f = fopen("energy.txt","a");
		double e = energy();
		fprintf(f,"%e %e %e\n",t, fabs((e-e_init)/e_init), tools_megno());
		fclose(f);
		printf("  Y = %.3f",tools_megno());
	}
}

tarch::la::Vector<NUMBER_OF_EULER_UNKNOWNS, double> peanoclaw::solver::euler3d::Cell::getUnknowns() const {
  tarch::la::Vector<NUMBER_OF_EULER_UNKNOWNS, double> unknowns;
  unknowns[0] = density();
  unknowns[1] = _data[1];
  unknowns[2] = _data[2];
  unknowns[3] = _data[3];
  unknowns[4] = energy();
  unknowns[5] = marker();

  return unknowns;
}

	int superMi::tuetee_mapping(Solid* mesh, double deltaE, double deltaT)
	{
		double oldEnergy = 0;
		double newEnergy = 0;
		int error = 0;
		bool flag = false;

		error = energy(mesh, &oldEnergy, TUETEE);
		std::cout << "Initial Tuette energy: " << oldEnergy << std::endl;
		for (int i = 0; (i < 100000) && !flag; i++){
			error = gradient(mesh, TUETEE);
			error = absolute_derivative(mesh);
			error = update_mesh(mesh, deltaT);
			error = energy(mesh, &newEnergy, TUETEE);
			if (i % 1000 == 0) std::cout << i + 1 << ": " << newEnergy << std::endl;
			error = check_energy_change(&oldEnergy, &newEnergy, deltaE, &flag);
		}

		return 0;
	}

void caller()
{
    int i;
    float *d;
    InteractionMatrix *jMatrices;
    d = (float*)malloc(sizeof(float)*3*3*15);
    for (i = 0; i < 15*9; i++) d[i] = i;
    jMatrices = (InteractionMatrix *)(d);
    printf("d=%p,jMatrices=%p\n",d,jMatrices);
    energy(jMatrices[3]);
}

/*
   This program computes first order directional derivatives 
   for the helmholtz energy function */
int main(int argc, char *argv[]) {
    int nf, n, j, l;
    double result1, result2;
    double q, jd, r;
    double *x, *bv;

    fprintf(stdout,"HELM-DIFF-EXAM (ADOL-C Example)\n\n");
    fprintf(stdout," # of independents/10 =? \n ");
    scanf("%d",&nf);

    /*--------------------------------------------------------------------------*/
    n = 10 * nf;                                            /* Initilizations */
    x   = (double*) malloc(n*sizeof(double));
    bv  = (double*) malloc(n*sizeof(double));

    r = 1.0/n;
    for (j=0; j<n; j++) {
        jd     = j;
        bv[j]  = 0.02*(1.0+fabs(sin(jd)));
        x[j]   = r*sqrt(1.0+jd);
    }

    /*--------------------------------------------------------------------------*/
    result2 = energy(n,x,bv);                                    /* basepoint */
    fprintf(stdout,"%14.6E -- energy\n",result2);

    /*--------------------------------------------------------------------------*/
    for (l=0; l<n; l++)                            /* directional derivatives */
    { x[l]    = x[l]+delta;
        result1 = energy(n,x,bv);
        x[l]    = x[l]-delta;
        q       = (result1-result2)/delta;
        fprintf(stdout,"%3d: %14.6E,  \n",l,q);
    }
    fprintf(stdout,"%14.6E -- energy\n",result2);

    free((char*) bv);
    free((char*) x);

    return 0;
}

bool uhfsolve::convergenceCriteria(){
    //Evaluate convergence conditions
    bool condition = true;
    if(iterations>5000){
        condition = false;
    }
    if(abs(energyPrev-energy())<tolerance){

        condition = false;
    }
    return condition;
}

//---------------------------------------------
double CpxCrvletPrtd::globalenergy()
{
  double lclsum = 0;
  vector< vector<int> >& c = _nx;
  for(int s=0; s<c.size(); s++)
	 for(int w=0; w<c[s].size(); w++)
		if(_owners[s][w]==mpirank())
		  lclsum += energy(_blocks[s][w]);
  double glbsum = 0;
  iC( MPI_Reduce((void*)(&lclsum), (void*)(&glbsum), 1, MPI_DOUBLE, MPI_SUM, 0, PETSC_COMM_WORLD) );
  return glbsum;
}

void problem_init(int argc, char* argv[]) {
    // Setup constants
    integrator	= WHFAST;
    dt 		= 0.001*2.*M_PI;			// initial timestep (in days)
    init_boxwidth(200);

    // Initial conditions

    {
        struct particle p = {.m=1.,.x=0,.y=0.,.z=0.,.vx=0,.vy=0.,.vz=0.};
        particles_add(p);
    }
    {
        double e = 0.999;
        struct particle p = {.m=0.,.x=0.01,.y=0.,.z=0.,.vx=0,.vy=0.*sqrt((1.+e)/(1.-e)),.vz=0.};
        particles_add(p);
    }
    tools_move_to_center_of_momentum();
    //problem_additional_forces 	= additional_forces;
    // Add megno particles
    //tools_megno_init(1e-16);  // N = 6 after this function call.
    system("rm -f *.txt");
    ei = energy();
}

void additional_forces() {
    particles[1].ax += 0.12/6.;
}

double energy() {
    double e_kin = 0.;
    double e_pot = 0.;
    struct particle pi = particles[1];
    e_kin += 0.5 * (pi.vx*pi.vx + pi.vy*pi.vy + pi.vz*pi.vz);
    struct particle pj = particles[0];
    double dx = pi.x - pj.x;
    double dy = pi.y - pj.y;
    double dz = pi.z - pj.z;
    e_pot -= G*pj.m/sqrt(dx*dx + dy*dy + dz*dz);
    return e_kin +e_pot;
}
int no =0;
void problem_output() {
    no++;
//	printf("%d\n", no);
    if (output_check(1000.*dt)) {
        output_timing();
    }
//	FILE* f = fopen("Y.txt","a+");
//	fprintf(f,"%e %e %e\n",t,(energy()-ei)/ei,tools_megno());
//	fclose(f);
}

CoordType SpatialModelMaximalRepulsion3D<CoordType>::computeBeta()
{
  const int n = 100;
  const int numPoints = this->getHardcoreDistances().getSize();
  const CoordType maxRadius = this->getTriMesh().equivalentRadius() / 50.0;
  RandomGenerator& randomGenerator = this->getRandomGenerator();
  Vertices<CoordType> vertices1;
  Vertices<CoordType> vertices2;
  CoordType sumDelta = 0.0;

  for (int i = 0; i < n; ++i)
  {
    vertices1 = SpatialModelHardcoreDistance3D<CoordType>::drawSample( numPoints );
    vertices2 = vertices1;
    vertices2.detach();
    int v = randomGenerator.uniformL( numPoints );
    moveVertex( vertices2, v, maxRadius );
    sumDelta += fabs( energy(vertices1) - energy(vertices2) );
  }

  return log( 20.0 ) / ( sumDelta/n );
}

TEST(podio, Basics) {
  auto store = podio::EventStore();
  // Adding
  auto& collection = store.create<ExampleHitCollection>("name");
  auto hit1 = collection.create(0.,0.,0.,0.); //initialize w/ value
  auto hit2 = collection.create(); //default initialize
  hit2.energy(12.5);
  // Retrieving
  const ExampleHitCollection* coll2(nullptr);
  bool success = store.get("name",coll2);
  const ExampleHitCollection* coll3(nullptr);
  if (store.get("wrongName",coll3) != false) success = false;
  EXPECT_EQ(true, success);
}

static PyObject *
energy_wrapped(PyObject *self, PyObject *args)
{
    float mass, e;

    if (!PyArg_ParseTuple(args, "f:energy", &mass))
    {
        return NULL;
    }
    
    e = energy(mass);
    
    return Py_BuildValue("f", e);
}

int main(int argc, char **argv){
    double* q_x = calloc(2, sizeof(double));
    double* q_y = calloc(2, sizeof(double));
    printf("Intial position in x : ");
    scanf("%lf",&q_x[0]);
    printf("Intial velocity in y : ");
    scanf("%lf",&q_y[1]);
    FILE *f = fopen("output.txt","w");
    FILE *f_1 = fopen("out3.txt","w");
    //q_x[0] = 1.0; // position
    q_x[1] = 0.0; // momentum 
    q_y[0] = 0.0;
    //q_y[1] = 1.0;

    int n = 10000;
    double h = (2*PI)/n;
    double int_E = energy(sqrt(pow(q_x[0],2) + pow(q_y[0],2)),q_x[1],q_y[1]);
    double int_L = angular_momentum(q_x[0],q_y[0],q_x[1],q_y[1]);
    fprintf(f,"x \t y \t R  \n");
    fprintf(f_1,"%s \t %s \t %s \t %s \t %s\n","Time","E","L","Deviation in E","Deviation in L");
    for (int i = 0;i<=n;i++){
        double radius = sqrt(pow(q_x[0],2) + pow(q_y[0],2));
        double theta = atan2(q_y[0],q_x[0]);
        double E = energy(radius,q_x[1],q_y[1]);
        double L = angular_momentum(q_x[0],q_y[0],q_x[1],q_y[1]);
        fprintf(f,"%f \t %f \t %f\n",q_x[0] ,q_y[0],radius);
        fprintf(f_1,"%.5e \t %.5e \t %.5e\t",i*h,E,L);
        fprintf(f_1,"%.6e \t %.16e\n",fabs(E-int_E), fabs(L-int_L));
        PERLF(q_x,force_x,h,radius);
        PERLF(q_y,force_y,h,radius);
        

    }
    fclose(f);
    return 0;

}