Merge pull request #87 from shishlo/master

The definition of Classical Radius of particle was fixed for heavy ions and electron
PyORBIT-Collaboration · Mar 5, 2024 · 8e54b39 · 8e54b39
2 parents dcceb7a + f85b4b7
commit 8e54b39
Show file tree

Hide file tree

Showing 17 changed files with 73 additions and 65 deletions.
diff --git a/Changes_to_PyORBIT2_Feb_March_2024.txt b/Changes_to_PyORBIT2_Feb_March_2024.txt
@@ -0,0 +1,17 @@
+==== Changes were made to the folowing files:
+src/orbit/Bunch.cc
+src/orbit/Bunch.hh
+src/orbit/wrap_bunch.cc
+src/spacecharge/BaseBoundary2D.cc
+src/spacecharge/ForceSolver2D.hh
+src/spacecharge/ForceSolverFFT2D.hh
+src/spacecharge/PoissonSolver2D.hh
+src/spacecharge/PoissonSolver3D.hh
+src/spacecharge/PoissonSolverFFT2D.hh
+src/spacecharge/PoissonSolverFFT3D.hh
+src/spacecharge/SpaceChargeCalc2p5D.cc
+src/spacecharge/SpaceChargeCalc2p5Drb.cc
+src/spacecharge/SpaceChargeCalc3D.cc
+src/spacecharge/SpaceChargeCalcSliceBySlice2D.cc
+src/spacecharge/SpaceChargeCalcUnifEllipse.cc
+src/spacecharge/SpaceChargeForceCalc2p5D.cc    
diff --git a/src/orbit/Bunch.cc b/src/orbit/Bunch.cc
@@ -157,6 +157,9 @@ void Bunch::setBunchAttribute(const std::string att_name, double att_val){
   if(att_name == "macro_size"){
     macroSizeForAll = att_val;
   }
+
+  //update classical radius of the particle
+  this->updateClassicalRadius();
 }
 
 void Bunch::setBunchAttribute(const std::string att_name, int att_val){
@@ -185,6 +188,9 @@ void Bunch::initBunchAttributes(const char* fileName){
   bunchAttr->doubleVal("charge",charge);
   bunchAttr->doubleVal("classical_radius",classicalRadius);
   bunchAttr->doubleVal("macro_size",macroSizeForAll);
+
+  //update classical radius of the particle
+  this->updateClassicalRadius();  
 
 	//read and set new ones
   std::vector<std::string> attr_names;
@@ -276,7 +282,10 @@ void Bunch::initBunchAttributes(const char* fileName){
   if(rank_MPI == 0){
     is.close();
   }
-
+
+  //update bunch attributes
+  this->updateClassicalRadius();
+
 	//set synchronous particle parameters from file
 	syncPart->readSyncPart(fileName);
 
@@ -767,19 +776,25 @@ double  Bunch::getMacroSize(){ return macroSizeForAll;}
 double Bunch::setMass(double val){
   mass = val;
   bunchAttr->doubleVal("mass",val);
+  this->updateClassicalRadius();
   return mass;
 }
 
 double Bunch::setCharge(double val){
   charge = val;
   bunchAttr->doubleVal("charge",val);
+  this->updateClassicalRadius();
   return charge;
 }
 
-double Bunch::setClassicalRadius(double val){
-  classicalRadius = val;
-  bunchAttr->doubleVal("classical_radius",val);
-  return classicalRadius;
+//Updates the classical radius of the particle according to mass and charge
+void Bunch::updateClassicalRadius(){
+	// charge in abs(electron charges)
+	//mass in GeV
+	double val = OrbitConst::classicalRadius_proton;
+	val *= OrbitConst::mass_proton/mass;
+	val *= charge*charge;
+	bunchAttr->doubleVal("classical_radius",val);
 }
 
 double  Bunch::setMacroSize(double val){

diff --git a/src/orbit/Bunch.hh b/src/orbit/Bunch.hh
@@ -113,7 +113,6 @@ public:
   double getMacroSize();
 
   double setMass(double mass);                // GeV
-  double setClassicalRadius(double clR);      // m
   double setCharge(double chrg);              // sign and value in abs(e-charge) only
   double setMacroSize(double mcrsz);
 
@@ -217,6 +216,9 @@ private:
   //array from ParticleAttributes class instance.
   //User is not supposed to use this method directly.
   double& getParticleAttributeVal(int ind, int attr_ind);
+
+  //Updates the classical radius of the particle according to mass and charge
+  void updateClassicalRadius();      // m   
 
 protected:
 

diff --git a/src/orbit/wrap_bunch.cc b/src/orbit/wrap_bunch.cc
@@ -498,39 +498,11 @@ namespace wrap_orbit_bunch{
     return Py_None;		
   }
 
-  //Sets or returns classicalRadius of the macro-particle in MeV
-  //  the action is depended on the number of arguments
-  //  classicalRadius() - returns classicalR
-  //  classicalRadius(value) - sets the new value
-  //this is implementation of the getClassicalRadius() and
-  //setClassicalRadius() methods of the Bunch class
+  //Returns classicalRadius of the particle in meters
   static PyObject* Bunch_classicalRadius(PyObject *self, PyObject *args){	
-		Bunch* cpp_bunch = (Bunch*) ((pyORBIT_Object *) self)->cpp_obj;		
-    //if nVars == 0 this is get classicalRadius
-    //if nVars == 1 this is set classicalRadius
-    int nVars = PyTuple_Size(args);
-
-    double val = 0.;
-
-    if(nVars == 0 ||  nVars == 1){
-      if(nVars == 0){
-        val = cpp_bunch->getClassicalRadius();
-      }
-      else{
-        //NO NEW OBJECT CREATED BY PyArg_ParseTuple! NO NEED OF Py_DECREF()
-        if(!PyArg_ParseTuple(	args,"d:classicalRadius",&val)){
-          error("PyBunch - classicalRadius(value) - value is needed");
-        }
-        cpp_bunch->setClassicalRadius(val);
-      }
-			return Py_BuildValue("d",val);
-    }
-    else{
-      error("PyBunch. You should call classicalRadius() or classicalRadius(value)");
-    }
-
-    Py_INCREF(Py_None);
-    return Py_None;		
+		Bunch* cpp_bunch = (Bunch*) ((pyORBIT_Object *) self)->cpp_obj;
+		double val = cpp_bunch->getClassicalRadius();
+		return Py_BuildValue("d",val);
   }
 
   //Sets or returns charge of the macro-particle in e-charge
@@ -1216,7 +1188,7 @@ namespace wrap_orbit_bunch{
     { "flag",                           Bunch_flag                          ,METH_VARARGS,"Returns flag(index) for particle with index"},
     { "ringwrap",                       Bunch_ringwrap                      ,METH_VARARGS,"Perform the ring wrap. Usage: ringwrap(ring_length)"},
     { "mass",                           Bunch_mass                          ,METH_VARARGS,"Set mass(value) or get mass() the mass of particle in MeV"},
-    { "classicalRadius",                Bunch_classicalRadius               ,METH_VARARGS,"Set and get a classical radius of particle in [m]"},
+    { "classicalRadius",                Bunch_classicalRadius               ,METH_VARARGS,"Returns a classical radius of particle in [m]"},
     { "charge",                         Bunch_charge                        ,METH_VARARGS,"Set charge(value) or get charge() the charge of particle in e-charge"},
     { "macroSize",                      Bunch_macroSize                     ,METH_VARARGS,"Set macroSize(value) or get macroSize() the charge of particle in e-charge"},
     { "initBunchAttr",                  Bunch_initBunchAttr                 ,METH_VARARGS,"Reads and initilizes bunch attributes from a bunch file"},

diff --git a/src/spacecharge/BaseBoundary2D.cc b/src/spacecharge/BaseBoundary2D.cc
@@ -1,4 +1,5 @@
 #include "BaseBoundary2D.hh"
+#include "OrbitConst.hh"
 
 #include <iostream>
 #include <cfloat>
@@ -9,8 +10,13 @@ const int BaseBoundary2D::IS_INSIDE    =  1;
 const int BaseBoundary2D::IS_OUTSIDE   = -1;
 const int BaseBoundary2D::TO_BE_KILLED =  0;
 
-const double BaseBoundary2D::PI = 3.14159265358979324;
+const double BaseBoundary2D::PI = OrbitConst::PI;
 
+/** 
+The BaseBoundary2D class defines a boundary geometry
+and calculates the potential created by charges on the boundary 
+surface.
+*/
 
 // Constructor
 BaseBoundary2D::BaseBoundary2D(int nPoints, int nModes): CppPyWrapper(NULL)

diff --git a/src/spacecharge/ForceSolver2D.hh b/src/spacecharge/ForceSolver2D.hh
@@ -19,9 +19,8 @@
 using namespace std;
 
 /** 
-  The ForceSolver2D class is used to define a boundary geometry
-  and to calculate the force created by charges on the boundary 
-	surface.
+  The ForceSolver2D class calculates 2D forces along X and Y axes
+  created by a 2D charge distribution.
 */
 
 class ForceSolver2D: public OrbitUtils::CppPyWrapper

diff --git a/src/spacecharge/ForceSolverFFT2D.hh b/src/spacecharge/ForceSolverFFT2D.hh
@@ -23,7 +23,8 @@
 using namespace std;
 
 /** 
-  The ForceSolverFFT2D class is used to calculate the force*/
+  The ForceSolverFFT2D class is used to calculate the force
+ */
 
 class ForceSolverFFT2D: public ForceSolver2D
 {

diff --git a/src/spacecharge/PoissonSolver2D.hh b/src/spacecharge/PoissonSolver2D.hh
@@ -19,9 +19,8 @@
 using namespace std;
 
 /** 
-  The PoissonSolver2D class is used to define a boundary geometry
-  and to calculate the potential created by charges on the boundary 
-	surface.
+  The PoissonSolver2D class calculates electrostatic
+  potential of a 2D charge distribution.
 */
 
 class PoissonSolver2D: public OrbitUtils::CppPyWrapper

diff --git a/src/spacecharge/PoissonSolver3D.hh b/src/spacecharge/PoissonSolver3D.hh
@@ -19,9 +19,8 @@
 using namespace std;
 
 /** 
-  The PoissonSolver3D class is used to define a boundary geometry
-  and to calculate the potential created by charges on the boundary 
-	surface.
+  The PoissonSolver3D class calculates 3D voltage distribution
+  created by 3D charge distribution.
 */
 
 class PoissonSolver3D: public OrbitUtils::CppPyWrapper

diff --git a/src/spacecharge/PoissonSolverFFT2D.hh b/src/spacecharge/PoissonSolverFFT2D.hh
@@ -22,9 +22,8 @@
 using namespace std;
 
 /** 
-  The PoissonSolverFFT2D class is used to define a boundary geometry
-  and to calculate the potential created by charges on the boundary 
-	surface.
+  The PoissonSolverFFT2D class calculates electrostatic
+  potential of a 2D charge distribution using 2D FFT approach.
 */
 
 class PoissonSolverFFT2D: public PoissonSolver2D

diff --git a/src/spacecharge/PoissonSolverFFT3D.hh b/src/spacecharge/PoissonSolverFFT3D.hh
@@ -23,9 +23,8 @@
 using namespace std;
 
 /** 
-  The PoissonSolverFFT3D class is used to define a boundary geometry
-  and to calculate the potential created by charges on the boundary 
-	surface.
+  The PoissonSolverFFT3D class calculates electrostatic
+  potential of a 3D charge distribution using 3D FFT approach.
 */
 
 class PoissonSolverFFT3D: public PoissonSolver3D

diff --git a/src/spacecharge/SpaceChargeCalc2p5D.cc b/src/spacecharge/SpaceChargeCalc2p5D.cc
@@ -94,7 +94,7 @@ void SpaceChargeCalc2p5D::trackBunch(Bunch* bunch, double length, BaseBoundary2D
 	}
 
 	SyncPart* syncPart = bunch->getSyncPart();	
-	double factor = 2*length*bunch->getClassicalRadius()*pow(bunch->getCharge(),2)/(pow(syncPart->getBeta(),2)*pow(syncPart->getGamma(),3));	
+	double factor = 2*length*bunch->getClassicalRadius()/(pow(syncPart->getBeta(),2)*pow(syncPart->getGamma(),3));	
 	//std::cout<<" debug totalMacrosize="<<totalMacrosize<<" factor="<<factor<<" z_step="<< z_step <<std::endl;	
 
 	factor = factor/(z_step*totalMacrosize);	

diff --git a/src/spacecharge/SpaceChargeCalc2p5Drb.cc b/src/spacecharge/SpaceChargeCalc2p5Drb.cc
@@ -125,7 +125,7 @@ void SpaceChargeCalc2p5Drb::trackBunch(Bunch* bunch, double length, double pipe_
 	poissonSolver->findPotential(rhoGrid,phiGrid);
 
 	SyncPart* syncPart = bunch->getSyncPart();	
-	double factor =  2*length*bunch->getClassicalRadius()*pow(bunch->getCharge(),2)/(pow(syncPart->getBeta(),2)*pow(syncPart->getGamma(),3));	
+	double factor =  2*length*bunch->getClassicalRadius()/(pow(syncPart->getBeta(),2)*pow(syncPart->getGamma(),3));	
 	//std::cout<<" debug totalMacrosize="<<totalMacrosize<<" factor="<<factor<<" z_step="<< z_step <<std::endl;	
 
 
@@ -138,7 +138,7 @@ void SpaceChargeCalc2p5Drb::trackBunch(Bunch* bunch, double length, double pipe_
 	double long_sc_factor_in = 1.0+2*log(pipe_radius/a_bunch);
 	double long_sc_factor_out = 2*log(pipe_radius);
 	double a_bunch_2 = a_bunch*a_bunch;
-	double long_sc_factor = - length*bunch->getClassicalRadius()*pow(bunch->getCharge(),2) * bunch->getMass()/(pow(syncPart->getGamma(),2));
+	double long_sc_factor = - length*bunch->getClassicalRadius()* bunch->getMass()/(pow(syncPart->getGamma(),2));
 	//std::cout<<" debug pipe_radius="<<pipe_radius<<" a_bunch="<<a_bunch<<std::endl;	
 	//std::cout<<" debug long_sc_factor_in="<<long_sc_factor_in<<" long_sc_factor_out="<<long_sc_factor_out<<std::endl;	
 	//std::cout<<" debug long_sc_factor="<<long_sc_factor<<std::endl;

diff --git a/src/spacecharge/SpaceChargeCalc3D.cc b/src/spacecharge/SpaceChargeCalc3D.cc
@@ -122,8 +122,8 @@ void SpaceChargeCalc3D::trackBunch(Bunch* bunch, double length){
 	//distance between bunches in the center of mass of the bunch
 	double lambda_cm = OrbitConst::c*beta*gamma/frequency_;
 
-	double trans_factor =  length*bunch->getClassicalRadius()*pow(bunch->getCharge(),2)/(pow(beta,2)*pow(gamma,2));	
-	double long_factor =  length*bunch->getClassicalRadius()*pow(bunch->getCharge(),2)*bunch->getMass();
+	double trans_factor =  length*bunch->getClassicalRadius()/(pow(beta,2)*pow(gamma,2));	
+	double long_factor =  length*bunch->getClassicalRadius()*bunch->getMass();
 
 	double x,y,z,ex,ey,ez;
 

diff --git a/src/spacecharge/SpaceChargeCalcSliceBySlice2D.cc b/src/spacecharge/SpaceChargeCalcSliceBySlice2D.cc
@@ -97,7 +97,7 @@ void SpaceChargeCalcSliceBySlice2D::trackBunch(Bunch* bunch, double length, Base
 	phiGrid3D->synchronizeMPI(bunch->getMPI_Comm_Local());
 
 	SyncPart* syncPart = bunch->getSyncPart();	
-	double factor = 2*length/rhoGrid3D->getStepZ()*bunch->getClassicalRadius()*pow(bunch->getCharge(),2)/(pow(syncPart->getBeta(),2)*pow(syncPart->getGamma(),3));	
+	double factor = 2*length/rhoGrid3D->getStepZ()*bunch->getClassicalRadius()/(pow(syncPart->getBeta(),2)*pow(syncPart->getGamma(),3));	
 
 	double x,y,z,ex,ey,ez;	
 

diff --git a/src/spacecharge/SpaceChargeCalcUnifEllipse.cc b/src/spacecharge/SpaceChargeCalcUnifEllipse.cc
@@ -76,8 +76,8 @@ void SpaceChargeCalcUnifEllipse::trackBunch(Bunch* bunch, double length){
 	//if there is nothing we give up
 	if(total_macrosize == 0.) return;
 
-	double trans_factor =  length*bunch->getClassicalRadius()*pow(bunch->getCharge(),2)/(pow(beta,2)*pow(gamma,2));	
-	double long_factor =  length*bunch->getClassicalRadius()*pow(bunch->getCharge(),2)*bunch->getMass();
+	double trans_factor =  length*bunch->getClassicalRadius()/(pow(beta,2)*pow(gamma,2));	
+	double long_factor =  length*bunch->getClassicalRadius()*bunch->getMass();
 
 	double x,y,z,ex,ey,ez;
 	for (int i = 0, n = bunch->getSize(); i < n; i++){

diff --git a/src/spacecharge/SpaceChargeForceCalc2p5D.cc b/src/spacecharge/SpaceChargeForceCalc2p5D.cc
@@ -89,7 +89,7 @@ void SpaceChargeForceCalc2p5D::trackBunch(Bunch* bunch, double length){
 	forceSolver->findForce(rhoGrid, forceGridX, forceGridY);
 
 	SyncPart* syncPart = bunch->getSyncPart();	
-	double factor = 2*length*bunch->getClassicalRadius()*pow(bunch->getCharge(),2)/(pow(syncPart->getBeta(),2)*pow(syncPart->getGamma(),3));	
+	double factor = 2*length*bunch->getClassicalRadius()/(pow(syncPart->getBeta(),2)*pow(syncPart->getGamma(),3));	
 
 	factor = factor/(z_step*totalMacrosize);	
 	double Lfactor = 0.;