diff --git a/ConeLFMode/ConeLFmode.cc b/ConeLFMode/ConeLFmode.cc index bacaddc..b7d38a7 100644 --- a/ConeLFMode/ConeLFmode.cc +++ b/ConeLFMode/ConeLFmode.cc @@ -7,8 +7,8 @@ // 04/21/2003 // // DESCRIPTION -// The base class for Python implementation of a field source. -// It should be sub-classed on Python level and implements +// The base class for Python implementation of a field source. +// It should be sub-classed on Python level and implements // getElectricField(x,y,z,t) and getMagneticField (x,y,z,t) methods. // The results of these methods will be available from the c++ level. // This is an example of embedding Python in C++ Orbit level. @@ -16,10 +16,10 @@ // /////////////////////////////////////////////////////////////////////////// #include "ConeLFmode.hh" -#include "MathPolynomial.hh" +#include -#include "orbit_mpi.hh" -#include "OrbitConst.hh" +#include +#include #include @@ -30,7 +30,7 @@ ConeLFmode::ConeLFmode(double P,double r_x,double r_y,double a_x,double a_y,doub { U=2*sqrt(P*OrbitConst::c*OrbitConst::permeability/OrbitConst::PI); - orient = new FieldOrientation(); + orient = new FieldOrientation(); Laser_lambda=la; k=2*OrbitConst::PI/Laser_lambda; rx=r_x; @@ -38,20 +38,20 @@ ConeLFmode::ConeLFmode(double P,double r_x,double r_y,double a_x,double a_y,doub ax=a_x; ay=a_y; - + //default orientation of laser field polarizations nEx=1;nEy=0;nEz=0; nHx=0;nHy=1;nHz=0; - - - + + + } ConeLFmode::~ConeLFmode() -{ +{ } @@ -65,26 +65,26 @@ ConeLFmode::~ConeLFmode() void ConeLFmode::setLaserFieldOrientation(double x0, double y0, double z0, double kx, double ky, double kz, double mx, double my, double mz, - double n_Ex, double n_Ey, double n_Ez) + double n_Ex, double n_Ey, double n_Ez) { double a=sqrt(n_Ex*n_Ex+n_Ey*n_Ey+n_Ez*n_Ez); orient->setCoefficients(x0, y0, z0,kx, ky, kz, mx, my, mz); - - nEx=n_Ex/a; - nEy=n_Ey/a; + + nEx=n_Ex/a; + nEy=n_Ey/a; nEz=n_Ez/a; //Vector of polarization of electric field - + double n_Hx=nEz*ky-nEy*kz; double n_Hy=nEx*kz-nEz*kx; double n_Hz=nEy*kx-nEx*ky; - + double b=sqrt(n_Hx*n_Hx+n_Hy*n_Hy+n_Hz*n_Hz); - + nHx=n_Hx/b; nHy=n_Hy/b; nHz=n_Hz/b; - + } @@ -97,10 +97,10 @@ double b=sqrt(n_Hx*n_Hx+n_Hy*n_Hy+n_Hz*n_Hz); double ConeLFmode::HermiteGaussianOmega(double m,double x, double y, double z, double px, double py, double pz, double t){ - - + + double a=OrbitConst::c/sqrt(m*m+px*px+py*py+pz*pz); - + double vx=px*a; double vy=py*a; double vz=pz*a; @@ -121,10 +121,10 @@ double ConeLFmode::HermiteGaussianOmega(double m,double x, double y, double z, d -void ConeLFmode::getLaserElectricMagneticField(double x, double y, double z, double t, +void ConeLFmode::getLaserElectricMagneticField(double x, double y, double z, double t, tcomplex& E_x, tcomplex& E_y, tcomplex& E_z, tcomplex& H_x, tcomplex& H_y, tcomplex& H_z){ - + orient->OrientCoordinates(x,y,z); @@ -137,7 +137,7 @@ tcomplex H=E/OrbitConst::c; E_x=E*nEx; E_y=E*nEy; E_z=E*nEz; H_x=H*nHx; H_y=H*nHy; H_z=H*nHz; - + } @@ -148,10 +148,10 @@ tcomplex H=E/OrbitConst::c; double ConeLFmode::getFrequencyOmega(double m, double x, double y, double z, double px, double py, double pz, double t){ - + orient->OrientCoordinates(x,y,z); orient->OrientVector2(px,py,pz); - + return HermiteGaussianOmega(m,x,y,z,px,py,pz,t); @@ -169,19 +169,19 @@ double ConeLFmode::getFrequencyOmega(double m, double x, double y, double z, dou tcomplex ConeLFmode::getNonOrientedU(double x, double y, double z, double t){ - - + + double rxax= 1/(rx+ax*z); double ryay= 1/(ry+ay*z); - + tcomplex p=tcomplex(0.,k*(z-t*OrbitConst::c+x*x*ax*rxax/2+y*y*ay*ryay/2)); - + return sqrt(rxax*ryay)*exp(-pow(x*rxax,2)-pow(y*ryay,2)-p); - - + + } diff --git a/ConeLFMode/wrap_cone_lf_mode.cc b/ConeLFMode/wrap_cone_lf_mode.cc index 29b7ba7..7ffd6ba 100644 --- a/ConeLFMode/wrap_cone_lf_mode.cc +++ b/ConeLFMode/wrap_cone_lf_mode.cc @@ -1,7 +1,6 @@ -#include "orbit_mpi.hh" -#include "pyORBIT_Object.hh" +#include +#include -#include "wrap_utils.hh" #include "wrap_cone_lf_mode.hh" #include @@ -12,7 +11,6 @@ #include "ConeLFmode.hh" using namespace OrbitUtils; -using namespace wrap_orbit_utils; namespace wrap_cone_lf_mode{ @@ -39,32 +37,32 @@ extern "C" { //initializator for python CppBaseFieldSource class //this is implementation of the __init__ method static int ConeLFmode_init(pyORBIT_Object *self, PyObject *args, PyObject *kwds){ - + double P=0; double rx=0; double ry=0; double a_x=0; double a_y=0; double lambda=0; - + if(!PyArg_ParseTuple( args,"dddddd:",&P,&rx,&ry,&a_x,&a_y,&lambda)){ error("LaserExternalEfects(P,rx,ry,ax,ay,lambda) - params. are needed"); } else { - + self->cpp_obj = new ConeLFmode(P,rx,ry,a_x,a_y,lambda); ((ConeLFmode*) self->cpp_obj)->setPyWrapper((PyObject*) self); } return 0; } - - - - - - + + + + + + static PyObject* ConeLFmode_setLaserFieldOrientation(PyObject *self, PyObject *args){ - + ConeLFmode* LaserField = (ConeLFmode*)((pyORBIT_Object*) self)->cpp_obj; double x0; double y0; @@ -82,23 +80,23 @@ extern "C" { if(!PyArg_ParseTuple( args,"dddddddddddd:",&x0, &y0, &z0, &kx, &ky, &kz, &mx, &my, &mz, &nEx, &nEy, &nEz)) {error("LaserExternalEfects - setLaserHalfAngle(comments is not written yet) - param. is needed");flag=false;} - + if(fabs(kx*mx+ky*my+kz*mz)/sqrt(kx*kx+ky*ky+kz*kz)/sqrt(mx*mx+my*my+mz*mz)>1.e-15) {error("Please be shure that kx*mx+ky*my+kz*mz==0");flag=false;} if(fabs(kx*nEx+ky*nEy+kz*nEz)/sqrt(kx*kx+ky*ky+kz*kz)/sqrt(nEx*nEx+nEy*nEy+nEz*nEz)>1.e-15) {error("Please be shure that kx*Ex+ky*Ey+kz*Ez==0");flag=false;} - - if(flag) + + if(flag) LaserField->setLaserFieldOrientation(x0,y0, z0, kx, ky, kz, mx, my, mz, nEx, nEy, nEz); - - + + Py_INCREF(Py_None); - return Py_None; - } - - - + return Py_None; + } + + + /* void LasStripExternalEffects::setLaserHalfAngle(double a); double LasStripExternalEffects::getLaserHalfAngle(); @@ -107,59 +105,59 @@ extern "C" { void LasStripExternalEffects::setLaser_lambda(double a); double LasStripExternalEffects::getLaser_lambda(); - + static PyObject* ConeLFmode_setLaserHalfAngle(PyObject *self, PyObject *args){ ConeLFmode* LaserField = (ConeLFmode*)((pyORBIT_Object*) self)->cpp_obj; - double half_Angle; + double half_Angle; if(!PyArg_ParseTuple( args,"d:",&half_Angle)){ error("LaserExternalEfects - setLaserHalfAngle(half_Angle) - param. is needed"); } else LaserField->setLaserHalfAngle(half_Angle); - - + + Py_INCREF(Py_None); - return Py_None; - } + return Py_None; + } static PyObject* ConeLFmode_setLaserPower(PyObject *self, PyObject *args){ ConeLFmode* LaserField = (ConeLFmode*)((pyORBIT_Object*) self)->cpp_obj; - double LaserPower; + double LaserPower; if(!PyArg_ParseTuple( args,"d:",&LaserPower)){ error("LaserExternalEfects - setLaserPower(LaserPower) - param. is needed"); } else LaserField->setLaserPower(LaserPower); - - + + Py_INCREF(Py_None); - return Py_None; + return Py_None; } static PyObject* ConeLFmode_setLaser_lambda(PyObject *self, PyObject *args){ ConeLFmode* LaserField = (ConeLFmode*)((pyORBIT_Object*) self)->cpp_obj; - double Laser_lambda; + double Laser_lambda; if(!PyArg_ParseTuple( args,"d:",&Laser_lambda)){ error("LaserExternalEfects - setLaser_lambda(Laser_lambda) - param. is needed"); } else LaserField->setLaserPower(Laser_lambda); - - + + Py_INCREF(Py_None); - return Py_None; + return Py_None; } - - - + + + static PyObject* ConeLFmode_setLaserFieldOrientation(PyObject *self, PyObject *args){ ConeLFmode* LaserField = (ConeLFmode*)((pyORBIT_Object*) self)->cpp_obj; - double Laser_lambda; + double Laser_lambda; if(!PyArg_ParseTuple( args,"d:",&Laser_lambda)){ error("LaserExternalEfects - setLaser_lambda(Laser_lambda) - param. is needed"); } else LaserField->setLaserPower(Laser_lambda); - - + + Py_INCREF(Py_None); - return Py_None; + return Py_None; } @@ -167,28 +165,28 @@ extern "C" { static PyObject* ConeLFmode_getLaserHalfAngle(PyObject *self, PyObject *args){ ConeLFmode* LaserField = (ConeLFmode*)((pyORBIT_Object*) self)->cpp_obj; - - return Py_BuildValue("d",LaserField->getLaserHalfAngle()); - } - + + return Py_BuildValue("d",LaserField->getLaserHalfAngle()); + } + static PyObject* ConeLFmode_getLaserPower(PyObject *self, PyObject *args){ ConeLFmode* LaserField = (ConeLFmode*)((pyORBIT_Object*) self)->cpp_obj; - - return Py_BuildValue("d",LaserField->getLaserPower()); - } + + return Py_BuildValue("d",LaserField->getLaserPower()); + } static PyObject* ConeLFmode_getLaser_lambda(PyObject *self, PyObject *args){ ConeLFmode* LaserField = (ConeLFmode*)((pyORBIT_Object*) self)->cpp_obj; - - return Py_BuildValue("d",LaserField->getLaser_omega()); - } + + return Py_BuildValue("d",LaserField->getLaser_omega()); + } */ - - - + + + //----------------------------------------------------- //destructor for python PyBaseFieldSource class (__del__ method). @@ -196,7 +194,7 @@ extern "C" { static void ConeLFmode_del(pyORBIT_Object* self){ //std::cerr<<"The CppBaseFieldSource __del__ has been called!"<cpp_obj); - self->ob_type->tp_free((PyObject*)self); + Py_TYPE(self)->tp_free((PyObject*)self); } // defenition of the methods of the python PyBaseFieldSource wrapper class @@ -222,8 +220,7 @@ extern "C" { //new python PyBaseFieldSource wrapper type definition static PyTypeObject pyORBIT_ConeLFmode_Type = { - PyObject_HEAD_INIT(NULL) - 0, /*ob_size*/ + PyVarObject_HEAD_INIT(NULL, 0) "ConeLFmode", /*tp_name*/ sizeof(pyORBIT_Object), /*tp_basicsize*/ 0, /*tp_itemsize*/ @@ -261,7 +258,7 @@ extern "C" { (initproc) ConeLFmode_init, /* tp_init */ 0, /* tp_alloc */ ConeLFmode_new, /* tp_new */ - }; + }; //-------------------------------------------------- //Initialization function of the pyPyBaseFieldSource class diff --git a/ContinuumSS/ContinuumSS.cc b/ContinuumSS/ContinuumSS.cc index 4f30ade..bbef411 100644 --- a/ContinuumSS/ContinuumSS.cc +++ b/ContinuumSS/ContinuumSS.cc @@ -7,26 +7,26 @@ // 06/27/2008 // // DESCRIPTION -// The base class for C++ implementation of a external effects -// during the transport of particles through the external field. +// The base class for C++ implementation of a external effects +// during the transport of particles through the external field. // It should be sub-classed on Python level and implement // setupEffects(Bunch* bunch) // finalizeEffects(Bunch* bunch) -// applyEffects(Bunch* bunch, int index, -// double* y_in_vct, double* y_out_vct, -// double t, double t_step, +// applyEffects(Bunch* bunch, int index, +// double* y_in_vct, double* y_out_vct, +// double t, double t_step, // OrbitUtils::BaseFieldSource* fieldSource) // methods. // The results of these methods will be available from the c++ level. // This is an example of embedding Python in C++ Orbit level. // -// -// -// -// +// +// +// +// // /////////////////////////////////////////////////////////////////////////// -#include "orbit_mpi.hh" +#include #include #include #include @@ -35,9 +35,9 @@ #include "ContinuumSS.hh" -#include "RungeKuttaTracker.hh" -#include "OrbitConst.hh" -#include "LorentzTransformationEM.hh" +#include +#include +//#include "LorentzTransformationEM.hh" #include "FieldRotation.hh" @@ -56,16 +56,16 @@ using namespace OrbitUtils; ContinuumSS::ContinuumSS(BaseLaserFieldSource* BaseLaserField,TDMcontinuum* cont) { - + setName("unnamed"); - + LaserField = BaseLaserField; - levels = cont->get_levels(); + levels = cont->get_levels(); + - dE = new double[levels + 1]; TDM = new double[levels + 1]; @@ -77,15 +77,15 @@ ContinuumSS::ContinuumSS(BaseLaserFieldSource* BaseLaserField,TDMcontinuum* cont delE = fabs(dE[2] - dE[1]); - - + + //allocating memory for koefficients of 4-th order Runge-Kutta method and other koeeficients of the master equation k_RungeKutt=new tcomplex*[levels+1]; for (int i=0;igetPyWrapper() == NULL){ delete LaserField; - } else { + } else { Py_XDECREF(LaserField->getPyWrapper()); } @@ -132,20 +132,20 @@ ContinuumSS::~ContinuumSS() void ContinuumSS::CalcPopulations(int i, Bunch* bunch) { - + for (int m=0; mattArr(i)[m] = pow(abs(a(i,m)),2); } - + /* double sum = 0; for (int m=1; mgetSize()]; for(int i=0; igetSize();i++) for(int n=1; nhasParticleAttributes("Amplitudes")==0) { std::map part_attr_dict; part_attr_dict["size"] = 2*levels+2; bunch->addParticleAttributes("Amplitudes",part_attr_dict); - + for (int i=0; igetSize();i++) bunch->getParticleAttributes("Amplitudes")->attValue(i,0) = 1; } - + if(bunch->hasParticleAttributes("Populations")==0) { std::map part_attr_dict; part_attr_dict["size"] = levels+1; bunch->addParticleAttributes("Populations",part_attr_dict); } - - + + AmplAttr = bunch->getParticleAttributes("Amplitudes"); PopAttr = bunch->getParticleAttributes("Populations"); - - - + + + for (int i=0; igetSize();i++) CalcPopulations(i, bunch); } - - + + void ContinuumSS::finalizeEffects(Bunch* bunch){ - + for (int i=0; igetSize();i++) CalcPopulations(i, bunch); - - for (int i=0;iregion(x,y,z)) { GetParticleFrameFields(i, t, t_step,bunch,fieldSource); - - // This function gives parameters Ez_stat Ex_las[1...3] Ey_las[1...3] Ez_las[1...3] t_part omega_part part_t_step - //in atomic units in frame of particle - GetParticleFrameParameters(i,t,t_step,bunch); + + // This function gives parameters Ez_stat Ex_las[1...3] Ey_las[1...3] Ez_las[1...3] t_part omega_part part_t_step + //in atomic units in frame of particle + GetParticleFrameParameters(i,t,t_step,bunch); // This function provides step solution for density matrix using rk4 method - AmplSolver4step(i,bunch); + AmplSolver4step(i,bunch); + - -// CalcPopulations(i, bunch); +// CalcPopulations(i, bunch); PopAttr->attArr(i)[0] = pow(abs(a(i,0)),2); @@ -261,38 +261,38 @@ if(LaserField->region(x,y,z)) { void ContinuumSS::GetParticleFrameFields(int i,double t, double t_step, Bunch* bunch, BaseFieldSource* fieldSource) { - + double E_x,E_y,E_z,absE; - + for (int j=0; j<3;j++) { - + exp_phasa[j] = LaserField->getLaserEMField(x0+j*(x-x0)/2,y0+j*(y-y0)/2,z0+j*(z-z0)/2, t+j*t_step/2,Ex_las[j],Ey_las[j],Ez_las[j],Bx_las[j],By_las[j],Bz_las[j]); - - - LorentzTransformationEM::complex_transform(bunch->getMass(), - px0,py0,pz0, - Ex_las[j],Ey_las[j],Ez_las[j], - Bx_las[j],By_las[j],Bz_las[j]); - - + + //LorentzTransformationEM::complex_transform(bunch->getMass(), + // px0,py0,pz0, + // Ex_las[j],Ey_las[j],Ez_las[j], + // Bx_las[j],By_las[j],Bz_las[j]); + + + if(j==0) absE = sqrt(pow(abs(Ex_las[j]),2)+pow(abs(Ey_las[j]),2)+pow(abs(Ez_las[j]),2)); - - - E_las[j] = absE*exp_phasa[j]; - + + E_las[j] = absE*exp_phasa[j]; + + } - - - + + + } @@ -303,8 +303,8 @@ void ContinuumSS::GetParticleFrameFields(int i,double t, double t_step, Bunch* void ContinuumSS::GetParticleFrameParameters(int i, double t,double t_step, Bunch* bunch) { - - + + double ta=2.418884326505e-17; //atomic unit of time double Ea=5.14220642e011; //Atomic unit of electric field double m=bunch->getMass(); @@ -316,8 +316,8 @@ double coeff=1./(gamma*ta); part_t_step=t_step*coeff; //time step in frame of particle (in atomic units) -for(int j=0;j<3;j++) - E_las[j]/=Ea; +for(int j=0;j<3;j++) + E_las[j]/=Ea; //cout.precision(20); //omega = gamma*ta*LaserField->getFrequencyOmega(m,x0,y0,z0,px0,py0,pz0,t); // frequensy of laser in particle frame (in atomic units) @@ -332,31 +332,31 @@ for(int j=0;j<3;j++) void ContinuumSS::AmplSolver4step(int i, Bunch* bunch) { - + tcomplex z,z1,z2,integral; double dt,JdEdt; - + for(int n=1; n #include "BaseLaserFieldSource.hh" #include "Stark.hh" -#include "ParticleAttributes.hh" +#include #include "TDMcontinuum.hh" diff --git a/ContinuumSS/TDMcontinuum.cc b/ContinuumSS/TDMcontinuum.cc index ae55ce6..4dfb473 100644 --- a/ContinuumSS/TDMcontinuum.cc +++ b/ContinuumSS/TDMcontinuum.cc @@ -7,19 +7,19 @@ // 04/21/2003 // // DESCRIPTION -// The base class for Python implementation of a field source. -// It should be sub-classed on Python level and implements +// The base class for Python implementation of a field source. +// It should be sub-classed on Python level and implements // getElectricField(x,y,z,t) and getMagneticField (x,y,z,t) methods. // The results of these methods will be available from the c++ level. // This is an example of embedding Python in C++ Orbit level. // // /////////////////////////////////////////////////////////////////////////// -#include "orbit_mpi.hh" -#include "BufferStore.hh" +#include +#include #include "TDMcontinuum.hh" -#include "LorentzTransformationEM.hh" +//#include "LorentzTransformationEM.hh" #include #include #include @@ -46,11 +46,11 @@ using namespace std; TDMcontinuum::TDMcontinuum(std::string addr_trans) -{ +{ int rank_MPI,size_MPI; ORBIT_MPI_Comm_size(MPI_COMM_WORLD, &size_MPI); - ORBIT_MPI_Comm_rank(MPI_COMM_WORLD, &rank_MPI); + ORBIT_MPI_Comm_rank(MPI_COMM_WORLD, &rank_MPI); ifstream file; ifstream files[MAX_NUMBER_OF_FILES]; double FF[MAX_NUMBER_OF_FILES]; @@ -62,102 +62,102 @@ TDMcontinuum::TDMcontinuum(std::string addr_trans) int k; - - + + if(rank_MPI == 0) { - + if((dp = opendir(addr_trans.c_str())) == NULL) cout << "Error(" << errno << ") opening " << addr_trans << endl; - + num_files = 0; levels=0; - while ((dirp = readdir(dp)) != NULL) + while ((dirp = readdir(dp)) != NULL) if(string(dirp->d_name)[0]==string("T")[0]) { num_files++; - + snprintf(addr_file,MAX_LENGTH_ADDRESS,"%s%s",addr_trans.c_str(),string(dirp->d_name).c_str()); files[num_files].open(addr_file); - if(num_files==1) { - while(!files[num_files].eof()) + if(num_files==1) { + while(!files[num_files].eof()) {if (levels==0) files[num_files]>>temp1>>temp; else files[num_files]>>temp2>>temp; levels++;} files[num_files].clear();files[num_files].close(); files[num_files].open(addr_file); } - + deltaE = temp2 - temp1; - + temp_str = string(dirp->d_name); temp_str.erase(0,6); temp_str.erase(temp_str.size()-4,4); FF[num_files] = atof(temp_str.c_str()); - + if(num_files == 1) {Fmin = FF[num_files]; Fmax = FF[num_files];} - else + else { if(FF[num_files] < Fmin) Fmin = FF[num_files]; if(FF[num_files] > Fmax) Fmax = FF[num_files]; } - + } - + closedir(dp); } - - - + + + ORBIT_MPI_Bcast(&levels,1 ,MPI_INT,0,MPI_COMM_WORLD); ORBIT_MPI_Bcast(&num_files,1 ,MPI_INT,0,MPI_COMM_WORLD); ORBIT_MPI_Bcast(&Fmin,1 ,MPI_DOUBLE,0,MPI_COMM_WORLD); ORBIT_MPI_Bcast(&Fmax,1 ,MPI_DOUBLE,0,MPI_COMM_WORLD); ORBIT_MPI_Bcast(&deltaE,1 ,MPI_DOUBLE,0,MPI_COMM_WORLD); - - - + + + dF = (Fmax - Fmin)/(num_files - 1); - + F = new double[num_files + 1]; E_F_tdm = new double*[levels+1]; for (int i=0;i>temp>>E_F_tdm[n][i]; if (n==1) dEmin[i] = 0.5 + temp; } - + for (int n=1;n>temp>>E_F_tdm[n][i]; if (n==1) dEmin[i] = 0.5 + temp; } + F[i] = FF[k]; FF[k] = 2*Fmax; - files[k].clear(); - files[k].close(); - + files[k].clear(); + files[k].close(); + } - + } - + for (int i=1;iFmax||E +#include #include diff --git a/ContinuumSS/wrap_continuum_ss.cc b/ContinuumSS/wrap_continuum_ss.cc index 0770f7f..b1bac41 100644 --- a/ContinuumSS/wrap_continuum_ss.cc +++ b/ContinuumSS/wrap_continuum_ss.cc @@ -1,7 +1,6 @@ -#include "orbit_mpi.hh" -#include "pyORBIT_Object.hh" +#include +#include -#include "wrap_utils.hh" #include "wrap_continuum_ss.hh" #include @@ -40,7 +39,7 @@ extern "C" { //initializator for python PyExternalEffects class //this is implementation of the __init__ method static int ContinuumSS_init(pyORBIT_Object *self, PyObject *args, PyObject *kwds){ - + @@ -49,21 +48,21 @@ extern "C" { if(!PyArg_ParseTuple( args,"OO:",&pyBaseLaserField,&pycontinuum)){ error("ContinuumSS(LaserField,StarkEffect) - params. are needed"); - } + } else { BaseLaserFieldSource* lfs = (BaseLaserFieldSource*) ((pyORBIT_Object*) pyBaseLaserField)->cpp_obj; TDMcontinuum* con = (TDMcontinuum*) ((pyORBIT_Object*) pycontinuum)->cpp_obj; self->cpp_obj = new ContinuumSS(lfs, con); ((ContinuumSS*) self->cpp_obj)->setPyWrapper((PyObject*) self); } - + return 0; } - - - - // name([name]) - sets or returns the name of the External Effeects class + + + + // name([name]) - sets or returns the name of the External Effeects class static PyObject* ContinuumSS_name(PyObject *self, PyObject *args){ ContinuumSS* cpp_ContinuumSS = (ContinuumSS*) ((pyORBIT_Object*) self)->cpp_obj; const char* name = NULL; @@ -75,20 +74,20 @@ extern "C" { cpp_ContinuumSS->setName(name_str); } return Py_BuildValue("s",cpp_ContinuumSS->getName().c_str()); - } - - - - - - + } + + + + + + //----------------------------------------------------- //destructor for python PyExternalEffects class (__del__ method). //----------------------------------------------------- static void ContinuumSS_del(pyORBIT_Object* self){ //std::cerr<<"The DensityMatrix __del__ has been called!"<cpp_obj); - self->ob_type->tp_free((PyObject*)self); + Py_TYPE(self)->tp_free((PyObject*)self); } // defenition of the methods of the python PyExternalEffects wrapper class @@ -108,8 +107,7 @@ extern "C" { //new python PyExternalEffects wrapper type definition static PyTypeObject pyORBIT_ContinuumSS_Type = { - PyObject_HEAD_INIT(NULL) - 0, /*ob_size*/ + PyVarObject_HEAD_INIT(NULL, 0) "ContinuumSS", /*tp_name*/ sizeof(pyORBIT_Object), /*tp_basicsize*/ 0, /*tp_itemsize*/ @@ -147,10 +145,10 @@ extern "C" { (initproc) ContinuumSS_init, /* tp_init */ 0, /* tp_alloc */ ContinuumSS_new, /* tp_new */ - }; + }; + - //-------------------------------------------------- //Initialization function of the pyPyExternalEffects class //It will be called from Bunch wrapper initialization @@ -159,7 +157,7 @@ extern "C" { if (PyType_Ready(&pyORBIT_ContinuumSS_Type) < 0) return; Py_INCREF(&pyORBIT_ContinuumSS_Type); PyModule_AddObject(module, "ContinuumSS", (PyObject *)&pyORBIT_ContinuumSS_Type); - + } #ifdef __cplusplus diff --git a/ContinuumSS/wrap_tdm_continuum.cc b/ContinuumSS/wrap_tdm_continuum.cc index afad378..06f217f 100644 --- a/ContinuumSS/wrap_tdm_continuum.cc +++ b/ContinuumSS/wrap_tdm_continuum.cc @@ -1,7 +1,6 @@ -#include "orbit_mpi.hh" -#include "pyORBIT_Object.hh" +#include +#include -#include "wrap_utils.hh" #include "wrap_tdm_continuum.hh" #include @@ -12,7 +11,6 @@ #include "TDMcontinuum.hh" //using namespace OrbitUtils; -using namespace wrap_orbit_utils; namespace wrap_tdm_continuum{ @@ -45,28 +43,28 @@ extern "C" { if(!PyArg_ParseTuple( args,"s:",&addressEG)){ error("TDMcontinuum(address,n1,n2,m) - params. are needed"); - } + } else { - + std::string addr(addressEG); self->cpp_obj = new TDMcontinuum(addr); ((TDMcontinuum*) self->cpp_obj)->setPyWrapper((PyObject*) self); } - + return 0; - - + + } - - - - + + + + static PyObject* TDMcontinuum_setField_returndE(PyObject *self, PyObject *args){ TDMcontinuum* cpp_TDMcontinuum = (TDMcontinuum*)((pyORBIT_Object*) self)->cpp_obj; - - + + int nVars = PyTuple_Size(args); double val; double mass; @@ -79,7 +77,7 @@ extern "C" { double px; double py; double pz; - + //NO NEW OBJECT CREATED BY PyArg_ParseTuple! NO NEED OF Py_DECREF() if(!PyArg_ParseTuple( args,"dddddddddd:",&mass,&E_x,&E_y,&E_z,&B_x,&B_y,&B_z,&px,&py,&pz)) @@ -88,12 +86,12 @@ extern "C" { val=cpp_TDMcontinuum->setField_returndE(mass,E_x,E_y,E_z,B_x,B_y,B_z,px,py,pz); return Py_BuildValue("d",val); } - - - - - + + + + + //----------------------------------------------------- //destructor for python PyBaseFieldSource class (__del__ method). @@ -101,7 +99,7 @@ extern "C" { static void TDMcontinuum_del(pyORBIT_Object* self){ //std::cerr<<"The CppBaseFieldSource __del__ has been called!"<cpp_obj); - self->ob_type->tp_free((PyObject*)self); + Py_TYPE(self)->tp_free((PyObject*)self); } // defenition of the methods of the python PyBaseFieldSource wrapper class @@ -127,8 +125,7 @@ extern "C" { //new python PyBaseFieldSource wrapper type definition static PyTypeObject pyORBIT_TDMcontinuum_Type = { - PyObject_HEAD_INIT(NULL) - 0, /*ob_size*/ + PyVarObject_HEAD_INIT(NULL, 0) "TDMcontinuum", /*tp_name*/ sizeof(pyORBIT_Object), /*tp_basicsize*/ 0, /*tp_itemsize*/ @@ -166,7 +163,7 @@ extern "C" { (initproc) TDMcontinuum_init, /* tp_init */ 0, /* tp_alloc */ TDMcontinuum_new, /* tp_new */ - }; + }; //-------------------------------------------------- //Initialization function of the pyPyBaseFieldSource class diff --git a/DensityMatrix/DensityMatrix.cc b/DensityMatrix/DensityMatrix.cc index 8b6ae30..91bf960 100644 --- a/DensityMatrix/DensityMatrix.cc +++ b/DensityMatrix/DensityMatrix.cc @@ -7,26 +7,26 @@ // 06/27/2008 // // DESCRIPTION -// The base class for C++ implementation of a external effects -// during the transport of particles through the external field. +// The base class for C++ implementation of a external effects +// during the transport of particles through the external field. // It should be sub-classed on Python level and implement // setupEffects(Bunch* bunch) // finalizeEffects(Bunch* bunch) -// applyEffects(Bunch* bunch, int index, -// double* y_in_vct, double* y_out_vct, -// double t, double t_step, +// applyEffects(Bunch* bunch, int index, +// double* y_in_vct, double* y_out_vct, +// double t, double t_step, // OrbitUtils::BaseFieldSource* fieldSource) // methods. // The results of these methods will be available from the c++ level. // This is an example of embedding Python in C++ Orbit level. // -// -// -// -// +// +// +// +// // /////////////////////////////////////////////////////////////////////////// -#include "orbit_mpi.hh" +#include #include #include #include @@ -35,9 +35,9 @@ //#include #include "DensityMatrix.hh" -#include "RungeKuttaTracker.hh" -#include "OrbitConst.hh" -#include "LorentzTransformationEM.hh" +#include +#include +//#include "LorentzTransformationEM.hh" #include "FieldRotation.hh" @@ -75,7 +75,7 @@ DensityMatrix::DensityMatrix(BaseLaserFieldSource* BaseLaserField, Stark* Stark_ Parameter_resonance=par_res; int st = StarkEffect->getStates(); levels = StarkEffect->getStates()*(1+StarkEffect->getStates())*(1+2*StarkEffect->getStates())/6; - + //allocating memory for koefficients of 4-th order Runge-Kutta method and other koeeficients of the master equation @@ -95,20 +95,20 @@ E_i=new double[levels+1]; for(int n=1;ngetPyWrapper() == NULL){ delete LaserField; } else { Py_XDECREF(LaserField->getPyWrapper()); } - - + + if(StarkEffect->getPyWrapper() == NULL){ delete LaserField; } else { Py_XDECREF(StarkEffect->getPyWrapper()); } - + } void DensityMatrix::CalcPopulations(int i, Bunch* bunch) { - + PopAttr->attArr(i)[0] = 1; for (int j=1; jattArr(i)[j] = Re(i,j,j); PopAttr->attArr(i)[0] -= PopAttr->attArr(i)[j]; } - + } -void DensityMatrix::setupEffects(Bunch* bunch){ - - +void DensityMatrix::setupEffects(Bunch* bunch){ + + int_E=new tcomplex**[levels+1]; for (int i=0;igetSize()]; - + for(int i=0; igetSize();i++) for(int n=1; nhasParticleAttributes("Amplitudes")==0) { std::map part_attr_dict; part_attr_dict["size"] = 2*levels*levels+1; bunch->addParticleAttributes("Amplitudes",part_attr_dict); - + for (int i=0; igetSize();i++) bunch->getParticleAttributes("Amplitudes")->attValue(i,1) = 1; } - + if(bunch->hasParticleAttributes("Populations")==0) { std::map part_attr_dict; part_attr_dict["size"] = levels+1; bunch->addParticleAttributes("Populations",part_attr_dict); } - - - - + + + + AmplAttr = bunch->getParticleAttributes("Amplitudes"); PopAttr = bunch->getParticleAttributes("Populations"); for (int i=0; igetSize();i++) CalcPopulations(i, bunch); - + } - - + + void DensityMatrix::finalizeEffects(Bunch* bunch){ - + for (int i=0;ix(0)<<"\t"<y(0)<<"\t"<z(0)<<"\n"; - + } @@ -269,45 +269,45 @@ void DensityMatrix::applyEffectsForEach(Bunch* bunch, int i, void DensityMatrix::GetParticleFrameFields(int i,double t,double t_step, Bunch* bunch, BaseFieldSource* fieldSource) { - + double Ez; - - fieldSource->getElectricMagneticField(x0,y0,z0,t,Ex_stat,Ey_stat,Ez_stat,Bx_stat,By_stat,Bz_stat); - LorentzTransformationEM::transform(bunch->getMass(),px0,py0,pz0,Ex_stat,Ey_stat,Ez_stat,Bx_stat,By_stat,Bz_stat); - - if(LaserField->region(x,y,z)) { - + fieldSource->getElectricMagneticField(x0,y0,z0,t,Ex_stat,Ey_stat,Ez_stat,Bx_stat,By_stat,Bz_stat); + //LorentzTransformationEM::transform(bunch->getMass(),px0,py0,pz0,Ex_stat,Ey_stat,Ez_stat,Bx_stat,By_stat,Bz_stat); + + + if(LaserField->region(x,y,z)) { + for (int j=0; j<3;j++) { - + LaserField->getLaserEMField(x0+j*(x-x0)/2,y0+j*(y-y0)/2,z0+j*(z-z0)/2, t+j*t_step/2,Ex_las[j],Ey_las[j],Ez_las[j],Bx_las[j],By_las[j],Bz_las[j]); - - LorentzTransformationEM::complex_transform(bunch->getMass(), - px0,py0,pz0, - Ex_las[j],Ey_las[j],Ez_las[j], - Bx_las[j],By_las[j],Bz_las[j]); - - + + //LorentzTransformationEM::complex_transform(bunch->getMass(), + // px0,py0,pz0, + // Ex_las[j],Ey_las[j],Ez_las[j], + // Bx_las[j],By_las[j],Bz_las[j]); + + Ez = FieldRotation::RotateElectricFieldsV(Ex_stat,Ey_stat,Ez_stat,Ex_las[j],Ey_las[j],Ez_las[j]); - + } } - else - for(int j=0;j<3;j++) { + else + for(int j=0;j<3;j++) { Ex_las[j] = 0.; Ey_las[j] = 0.; Ez_las[j] = 0.; } - - Ex_stat=0; - Ey_stat=0; + + Ex_stat=0; + Ey_stat=0; Ez_stat=Ez; - - - + + + } @@ -318,8 +318,8 @@ void DensityMatrix::GetParticleFrameFields(int i,double t,double t_step, Bunch* void DensityMatrix::GetParticleFrameParameters(int i, double t,double t_step, Bunch* bunch) { - - + + double ta=2.418884326505e-17; //atomic unit of time double Ea=5.14220642e011; //Atomic unit of electric field double m=bunch->getMass(); @@ -333,14 +333,14 @@ part_t_step=t_step*coeff; //time step in frame of particle (in atomic units) -for(int j=0;j<3;j++) -{ - Ex_las[j]/=Ea; - Ey_las[j]/=Ea; +for(int j=0;j<3;j++) +{ + Ex_las[j]/=Ea; + Ey_las[j]/=Ea; Ez_las[j]/=Ea; } -Ex_stat/=Ea; -Ey_stat/=Ea; +Ex_stat/=Ea; +Ey_stat/=Ea; Ez_stat/=Ea; @@ -348,7 +348,7 @@ Ez_stat/=Ea; omega_part=gamma*ta*LaserField->getFrequencyOmega(m,x0,y0,z0,px0,py0,pz0,t); // frequensy of laser in particle frame (in atomic units) - + } @@ -365,50 +365,50 @@ omega_part=gamma*ta*LaserField->getFrequencyOmega(m,x0,y0,z0,px0,py0,pz0,t); // void DensityMatrix::AmplSolver4step(int i, Bunch* bunch) { - - - + + + tcomplex JdEdt, z,z1,z2,mu_x,mu_y,mu_z; double dt,sum, dE; - - - - StarkEffect->SetE(Ez_stat); - - - - - - + + + + StarkEffect->SetE(Ez_stat); + + + + + + for(int n=2; nfield_thresh[n]>Ez_stat && StarkEffect->field_thresh[m]>Ez_stat) { - + gamma_ij[n][m] = StarkEffect->getRelaxTransition(n,m,mu_x,mu_y,mu_z); dE = fabs(StarkEffect->En[n]-StarkEffect->En[m]); JdEdt = tcomplex(0.,dE*part_t_step); - + mu_Elas[n][m][1]=(mu_x*conj(Ex_las[1])+mu_y*conj(Ey_las[1])+mu_z*conj(Ez_las[1]))*(J/2.); cond[n][m]=fabs(dE - omega_part) < Parameter_resonance*abs(2.*mu_Elas[n][m][1]); - + if (cond[n][m]) { - - mu_Elas[n][m][0]=(mu_x*conj(Ex_las[0])+mu_y*conj(Ey_las[0])+mu_z*conj(Ez_las[0]))*(J/2.); + + mu_Elas[n][m][0]=(mu_x*conj(Ex_las[0])+mu_y*conj(Ey_las[0])+mu_z*conj(Ez_las[0]))*(J/2.); mu_Elas[n][m][2]=(mu_x*conj(Ex_las[2])+mu_y*conj(Ey_las[2])+mu_z*conj(Ez_las[2]))*(J/2.); - - z1=exp(int_E[n][m][i]); + + z1=exp(int_E[n][m][i]); z2=exp(JdEdt/2.); - - exp_mu_El[n][m][1]=z1*mu_Elas[n][m][0]; + + exp_mu_El[n][m][1]=z1*mu_Elas[n][m][0]; exp_mu_El[n][m][2]=z1*z2*mu_Elas[n][m][1]; exp_mu_El[n][m][3]=exp_mu_El[n][m][2]; exp_mu_El[n][m][4]=z1*z2*z2*mu_Elas[n][m][2]; - + } - + int_E[n][m][i] += JdEdt; - + } @@ -416,22 +416,22 @@ void DensityMatrix::AmplSolver4step(int i, Bunch* bunch) { - - + + for(int j=1; j<5; j++) for(int n=1; nfield_thresh[m]>Ez_stat && StarkEffect->field_thresh[n]>Ez_stat){ - + for(int m=1; mfield_thresh[m]>Ez_stat && StarkEffect->field_thresh[n]>Ez_stat){ + if (j==4) dt=part_t_step; else dt=part_t_step/2.; - + k_RungeKutt[n][m][j]=tcomplex(0); - for(int k=1;kfield_thresh[k]>Ez_stat) k_RungeKutt[n][m][j]+=exp_mu_El[n][k][j]*(dm(i,k,m)+k_RungeKutt[k][m][j-1]*dt); - for(int k=n+1;kfield_thresh[k]>Ez_stat) k_RungeKutt[n][m][j]-=conj(exp_mu_El[k][n][j])*(dm(i,k,m)+k_RungeKutt[k][m][j-1]*dt); - for(int k=1;kfield_thresh[k]>Ez_stat) k_RungeKutt[n][m][j]+=conj(exp_mu_El[m][k][j])*(dm(i,n,k)+k_RungeKutt[n][k][j-1]*dt); - for(int k=m+1;kfield_thresh[k]>Ez_stat) k_RungeKutt[n][m][j]-=exp_mu_El[k][m][j]*(dm(i,n,k)+k_RungeKutt[n][k][j-1]*dt); + for(int k=1;kfield_thresh[k]>Ez_stat) k_RungeKutt[n][m][j]+=exp_mu_El[n][k][j]*(dm(i,k,m)+k_RungeKutt[k][m][j-1]*dt); + for(int k=n+1;kfield_thresh[k]>Ez_stat) k_RungeKutt[n][m][j]-=conj(exp_mu_El[k][n][j])*(dm(i,k,m)+k_RungeKutt[k][m][j-1]*dt); + for(int k=1;kfield_thresh[k]>Ez_stat) k_RungeKutt[n][m][j]+=conj(exp_mu_El[m][k][j])*(dm(i,n,k)+k_RungeKutt[n][k][j-1]*dt); + for(int k=m+1;kfield_thresh[k]>Ez_stat) k_RungeKutt[n][m][j]-=exp_mu_El[k][m][j]*(dm(i,n,k)+k_RungeKutt[n][k][j-1]*dt); @@ -440,41 +440,41 @@ void DensityMatrix::AmplSolver4step(int i, Bunch* bunch) { if(n==m){ for(int k=m+1;kfield_thresh[k]>Ez_stat) k_RungeKutt[n][m][j]+=gamma_ij[k][m]*(dm(i,k,k)+k_RungeKutt[k][k][j-1]*dt); sum=0; for(int k=1;kfield_thresh[k]>Ez_stat) sum+=gamma_ij[m][k]; - + k_RungeKutt[n][m][j]-=(dm(i,n,m)+k_RungeKutt[n][m][j-1]*dt)*sum; } - + if(n!=m){ sum=0; for(int k=1;kfield_thresh[k]>Ez_stat) sum+=gamma_ij[m][k]; for(int k=1;kfield_thresh[k]>Ez_stat) sum+=gamma_ij[n][k]; k_RungeKutt[n][m][j]-=(dm(i,n,m)+k_RungeKutt[n][m][j-1]*dt)*sum/2.; } - - + + k_RungeKutt[n][m][j]-=(StarkEffect->Gamman[n]+StarkEffect->Gamman[m])*(dm(i,n,m)+k_RungeKutt[n][m][j-1]*dt)/2.; - + k_RungeKutt[m][n][j]=conj(k_RungeKutt[n][m][j]); - - + + } - - - - - + + + + + for(int n=1;nfield_thresh[m]>Ez_stat && StarkEffect->field_thresh[n]>Ez_stat) { - - z=(k_RungeKutt[n][m][1]+2.*k_RungeKutt[n][m][2]+2.*k_RungeKutt[n][m][3]+k_RungeKutt[n][m][4])/6.; + + z=(k_RungeKutt[n][m][1]+2.*k_RungeKutt[n][m][2]+2.*k_RungeKutt[n][m][3]+k_RungeKutt[n][m][4])/6.; Re(i,n,m)+=z.real()*part_t_step; Im(i,n,m)+=z.imag()*part_t_step; - + Re(i,m,n)=Re(i,n,m); Im(i,m,n)=-Im(i,n,m); } @@ -484,8 +484,8 @@ void DensityMatrix::AmplSolver4step(int i, Bunch* bunch) { Re(i,n,m) = 0; Im(i,n,m) = 0; } - + } } - + diff --git a/DensityMatrix/DensityMatrix.hh b/DensityMatrix/DensityMatrix.hh index b1501ff..0e73738 100644 --- a/DensityMatrix/DensityMatrix.hh +++ b/DensityMatrix/DensityMatrix.hh @@ -5,10 +5,10 @@ #include "Python.h" -#include "ExternalEffects.hh" +#include #include "BaseLaserFieldSource.hh" #include "Stark.hh" -#include "ParticleAttributes.hh" +#include using namespace TrackerRK4; diff --git a/DensityMatrix/wrap_density_matrix.cc b/DensityMatrix/wrap_density_matrix.cc index 130e254..79f2cd1 100644 --- a/DensityMatrix/wrap_density_matrix.cc +++ b/DensityMatrix/wrap_density_matrix.cc @@ -1,7 +1,6 @@ -#include "orbit_mpi.hh" -#include "pyORBIT_Object.hh" +#include +#include -#include "wrap_utils.hh" #include "wrap_density_matrix.hh" #include @@ -39,32 +38,32 @@ extern "C" { //initializator for python PyExternalEffects class //this is implementation of the __init__ method static int DensityMatrix_init(pyORBIT_Object *self, PyObject *args, PyObject *kwds){ - - + + double par_res=0; PyObject* pyBaseLaserField=NULL; PyObject* pyStarkEffect=NULL; - + if(!PyArg_ParseTuple( args,"OOd:",&pyBaseLaserField,&pyStarkEffect,&par_res)) - {error("DensityMatrix(LaserField,StarkEffect,resonanse_parameter) - params. are needed");} + {error("DensityMatrix(LaserField,StarkEffect,resonanse_parameter) - params. are needed");} else { BaseLaserFieldSource* lfs = (BaseLaserFieldSource*) ((pyORBIT_Object*) pyBaseLaserField)->cpp_obj; Stark* Stark_ef = (Stark*) ((pyORBIT_Object*) pyStarkEffect)->cpp_obj; self->cpp_obj = new DensityMatrix(lfs, Stark_ef, par_res); ((DensityMatrix*) self->cpp_obj)->setPyWrapper((PyObject*) self); } - - + + return 0; } - - - - // name([name]) - sets or returns the name of the External Effeects class + + + + // name([name]) - sets or returns the name of the External Effeects class static PyObject* DensityMatrix_name(PyObject *self, PyObject *args){ DensityMatrix* cpp_DensityMatrix = (DensityMatrix*) ((pyORBIT_Object*) self)->cpp_obj; const char* name = NULL; @@ -76,32 +75,32 @@ extern "C" { cpp_DensityMatrix->setName(name_str); } return Py_BuildValue("s",cpp_DensityMatrix->getName().c_str()); - } - - - - - - - - - - - - - - - - - - + } + + + + + + + + + + + + + + + + + + //----------------------------------------------------- //destructor for python PyExternalEffects class (__del__ method). //----------------------------------------------------- static void DensityMatrix_del(pyORBIT_Object* self){ //std::cerr<<"The DensityMatrix __del__ has been called!"<cpp_obj); - self->ob_type->tp_free((PyObject*)self); + Py_TYPE(self)->tp_free((PyObject*)self); } // defenition of the methods of the python PyExternalEffects wrapper class @@ -120,8 +119,7 @@ extern "C" { //new python PyExternalEffects wrapper type definition static PyTypeObject pyORBIT_DensityMatrix_Type = { - PyObject_HEAD_INIT(NULL) - 0, /*ob_size*/ + PyVarObject_HEAD_INIT(NULL, 0) "DensityMatrix", /*tp_name*/ sizeof(pyORBIT_Object), /*tp_basicsize*/ 0, /*tp_itemsize*/ @@ -159,10 +157,10 @@ extern "C" { (initproc) DensityMatrix_init, /* tp_init */ 0, /* tp_alloc */ DensityMatrix_new, /* tp_new */ - }; + }; + - //-------------------------------------------------- //Initialization function of the pyPyExternalEffects class //It will be called from Bunch wrapper initialization @@ -171,7 +169,7 @@ extern "C" { if (PyType_Ready(&pyORBIT_DensityMatrix_Type) < 0) return; Py_INCREF(&pyORBIT_DensityMatrix_Type); PyModule_AddObject(module, "DensityMatrix", (PyObject *)&pyORBIT_DensityMatrix_Type); - + } #ifdef __cplusplus diff --git a/DiagonalDensityMatrix/DM_noLaserField.cc b/DiagonalDensityMatrix/DM_noLaserField.cc index 7388537..7a30d7f 100644 --- a/DiagonalDensityMatrix/DM_noLaserField.cc +++ b/DiagonalDensityMatrix/DM_noLaserField.cc @@ -7,20 +7,20 @@ // 06/27/2008 // // DESCRIPTION -// -// +// +// // /////////////////////////////////////////////////////////////////////////// -#include "orbit_mpi.hh" +#include #include #include #include #include #include "DM_noLaserField.hh" -#include "RungeKuttaTracker.hh" -#include "OrbitConst.hh" -#include "LorentzTransformationEM.hh" +#include +#include +//#include "LorentzTransformationEM.hh" #define pop(i,n) PopAttr->attArr(i)[n] //i-part index, n,m-attr index @@ -38,41 +38,41 @@ DM_noLaserField::DM_noLaserField(Stark* Starkef, int method) index = method; StarkEffect=Starkef; - + int st = StarkEffect->getStates(); levels= st*(1+st)*(1+2*st)/6; - + //allocating memory for koefficients of 4-th order Runge-Kutta method and other koeeficients of the master equation - + k_RungeKutt=new double*[levels+1]; for (int i=0;igetPyWrapper() != NULL){ Py_INCREF(StarkEffect->getPyWrapper()); } - + } DM_noLaserField::~DM_noLaserField() @@ -80,7 +80,7 @@ DM_noLaserField::~DM_noLaserField() for (int i=0;igetPyWrapper() == NULL){ delete StarkEffect; } else { @@ -90,16 +90,16 @@ DM_noLaserField::~DM_noLaserField() void DM_noLaserField::CalcPopulations(int i, Bunch* bunch) { PopAttr->attArr(i)[0] = 1; - for (int j=1; jattArr(i)[0] -= PopAttr->attArr(i)[j]; - + } -void DM_noLaserField::setupEffects(Bunch* bunch){ - +void DM_noLaserField::setupEffects(Bunch* bunch){ + prob = new double[bunch->getSize()]; bunch->setCharge(+1); @@ -108,59 +108,59 @@ void DM_noLaserField::setupEffects(Bunch* bunch){ prob[i] = (double)rand()/(double)RAND_MAX; // std::cout<hasParticleAttributes("Populations")!=0) if(bunch->getParticleAttributes("Populations")->getAttSize() != levels+1) bunch->removeAllParticleAttributes(); - - + + if(bunch->hasParticleAttributes("Populations")==0) { std::map part_attr_dict; part_attr_dict["size"] = levels+1; bunch->addParticleAttributes("Populations",part_attr_dict); - + for (int i=0; igetSize();i++) bunch->getParticleAttributes("Populations")->attValue(i,1) = 1; } - - + + PopAttr = bunch->getParticleAttributes("Populations"); - + for (int i=0; igetSize();i++) CalcPopulations(i, bunch); } - + void DM_noLaserField::finalizeEffects(Bunch* bunch){ - + delete [] prob; - + } -void DM_noLaserField::applyEffectsForEach(Bunch* bunch, int i, - double* y_in_vct, double* y_out_vct, - double t, double t_step, +void DM_noLaserField::applyEffectsForEach(Bunch* bunch, int i, + double* y_in_vct, double* y_out_vct, + double t, double t_step, OrbitUtils::BaseFieldSource* fieldSource, RungeKuttaTracker* tracker) { if(bunch->flag(i) == 0) { - + double** partCoordArr = bunch->coordArr(); - + y_in_vct[0] = partCoordArr[i][0]; y_in_vct[1] = partCoordArr[i][2]; y_in_vct[2] = partCoordArr[i][4]; y_in_vct[3] = partCoordArr[i][1]; y_in_vct[4] = partCoordArr[i][3]; - y_in_vct[5] = partCoordArr[i][5]; - - double mass = bunch->getMass(); + y_in_vct[5] = partCoordArr[i][5]; + + double mass = bunch->getMass(); double vpt = mass*mass + y_in_vct[3]*y_in_vct[3] + y_in_vct[4]*y_in_vct[4] + y_in_vct[5]*y_in_vct[5]; - + vpt = 299792458*t_step/sqrt(vpt); y_out_vct[0] = y_in_vct[0] + vpt*y_in_vct[3]; y_out_vct[1] = y_in_vct[1] + vpt*y_in_vct[4]; @@ -168,16 +168,16 @@ void DM_noLaserField::applyEffectsForEach(Bunch* bunch, int i, y_out_vct[3] = y_in_vct[3]; y_out_vct[4] = y_in_vct[4]; y_out_vct[5] = y_in_vct[5]; - + partCoordArr[i][0] = y_out_vct[0]; partCoordArr[i][2] = y_out_vct[1]; partCoordArr[i][4] = y_out_vct[2]; partCoordArr[i][1] = y_out_vct[3]; partCoordArr[i][3] = y_out_vct[4]; - partCoordArr[i][5] = y_out_vct[5]; - + partCoordArr[i][5] = y_out_vct[5]; + } - + x0 = y_in_vct[0];y0 = y_in_vct[1];z0 = y_in_vct[2]; px0 = y_in_vct[3];py0 = y_in_vct[4];pz0 = y_in_vct[5]; @@ -187,17 +187,17 @@ void DM_noLaserField::applyEffectsForEach(Bunch* bunch, int i, - // This function gives parameters Ez_stat Ex_las[1...3] Ey_las[1...3] Ez_las[1...3] - //in natural unts (Volt per meter) in the frame of particle + // This function gives parameters Ez_stat Ex_las[1...3] Ey_las[1...3] Ez_las[1...3] + //in natural unts (Volt per meter) in the frame of particle GetParticleFrameFields(i, t,t_step, bunch,fieldSource); - - // This function gives parameters Ez_stat Ex_las[1...3] Ey_las[1...3] Ez_las[1...3] t_part omega_part part_t_step - //in atomic units in frame of particle - GetParticleFrameParameters(i,t,t_step,bunch); + + // This function gives parameters Ez_stat Ex_las[1...3] Ey_las[1...3] Ez_las[1...3] t_part omega_part part_t_step + //in atomic units in frame of particle + GetParticleFrameParameters(i,t,t_step,bunch); // This function provides step solution for density matrix using rk4 method - AmplSolver4step(t_step, i,bunch); - + AmplSolver4step(t_step, i,bunch); + CalcPopulations(i, bunch); //cout<x(0)<<"\t"<y(0)<<"\t"<z(0)<<"\t"<attArr(i)[0]<<"\t" <attArr(i)[12] + PopAttr->attArr(i)[13]+PopAttr->attArr(i)[7] + PopAttr->attArr(i)[8] + PopAttr->attArr(i)[1]<<"\n"; @@ -205,34 +205,34 @@ void DM_noLaserField::applyEffectsForEach(Bunch* bunch, int i, void DM_noLaserField::GetParticleFrameFields(int i,double t,double t_step, Bunch* bunch, BaseFieldSource* fieldSource) { double Ez; - - fieldSource->getElectricMagneticField(x0,y0,z0,t,Ex_stat,Ey_stat,Ez_stat,Bx_stat,By_stat,Bz_stat); - LorentzTransformationEM::transform(bunch->getMass(),px0,py0,pz0,Ex_stat,Ey_stat,Ez_stat,Bx_stat,By_stat,Bz_stat); - + + fieldSource->getElectricMagneticField(x0,y0,z0,t,Ex_stat,Ey_stat,Ez_stat,Bx_stat,By_stat,Bz_stat); + //LorentzTransformationEM::transform(bunch->getMass(),px0,py0,pz0,Ex_stat,Ey_stat,Ez_stat,Bx_stat,By_stat,Bz_stat); + Ez=sqrt(Ex_stat*Ex_stat+Ey_stat*Ey_stat+Ez_stat*Ez_stat); - - Ex_stat=0; - Ey_stat=0; + + Ex_stat=0; + Ey_stat=0; Ez_stat=Ez; } -void DM_noLaserField::GetParticleFrameParameters(int i, double t,double t_step, Bunch* bunch) { - +void DM_noLaserField::GetParticleFrameParameters(int i, double t,double t_step, Bunch* bunch) { + double ta=2.418884326505e-17; //atomic unit of time double Ea=5.14220642e011; //Atomic unit of electric field double m=bunch->getMass(); - + //This line calculates relyativistic factor-Gamma double gamma=sqrt(m*m+px0*px0+py0*py0+pz0*pz0)/m; double coeff=1./(gamma*ta); - + part_t_step=t_step*coeff; //time step in frame of particle (in atomic units) - t_part=t*coeff; - - - Ex_stat/=Ea; - Ey_stat/=Ea; - Ez_stat/=Ea; + t_part=t*coeff; + + + Ex_stat/=Ea; + Ey_stat/=Ea; + Ez_stat/=Ea; } @@ -240,12 +240,12 @@ void DM_noLaserField::GetParticleFrameParameters(int i, double t,double t_step, void DM_noLaserField::AmplSolver4step(double t_step, int i, Bunch* bunch) { double dt,sum; - - + + if(bunch->flag(i) == 0) { - - StarkEffect->SetE(Ez_stat); - + + StarkEffect->SetE(Ez_stat); + for(int n=2; nfield_thresh[n]>Ez_stat && StarkEffect->field_thresh[m]>Ez_stat) { @@ -254,51 +254,51 @@ void DM_noLaserField::AmplSolver4step(double t_step, int i, Bunch* bunch) { } } } - - + + for(int j=1; j<5; j++){ for(int m=1; mfield_thresh[m]>Ez_stat) { - + if (StarkEffect->field_thresh[m]>Ez_stat) { + if (j==4) dt=part_t_step; else dt=part_t_step/2.; - + k_RungeKutt[m][j]=0; - - for(int k=m+1;kfield_thresh[k]>Ez_stat) k_RungeKutt[m][j]+=gamma_ij[k][m]*(pop(i,k)+k_RungeKutt[k][j-1]*dt); + + for(int k=m+1;kfield_thresh[k]>Ez_stat) k_RungeKutt[m][j]+=gamma_ij[k][m]*(pop(i,k)+k_RungeKutt[k][j-1]*dt); sum=0;for(int k=1;kfield_thresh[k]>Ez_stat) sum+=gamma_ij[m][k]; - + sum+=StarkEffect->Gamman[m]; k_RungeKutt[m][j]-=(pop(i,m)+k_RungeKutt[m][j-1]*dt)*sum; - - + + } } } - + for(int m=1;mfield_thresh[m]>Ez_stat) { - pop(i,m) += part_t_step*(k_RungeKutt[m][1]+2.*k_RungeKutt[m][2]+2.*k_RungeKutt[m][3]+k_RungeKutt[m][4])/6.; + pop(i,m) += part_t_step*(k_RungeKutt[m][1]+2.*k_RungeKutt[m][2]+2.*k_RungeKutt[m][3]+k_RungeKutt[m][4])/6.; } else { pop(i,m) = 0; } - + } - + } - - - + + + if((index == 1)&&(pop(i,0) >= prob[i])) bunch->recoverParticle(i); - - + + } diff --git a/DiagonalDensityMatrix/DM_noLaserField.hh b/DiagonalDensityMatrix/DM_noLaserField.hh index d947af7..d8dd9ab 100644 --- a/DiagonalDensityMatrix/DM_noLaserField.hh +++ b/DiagonalDensityMatrix/DM_noLaserField.hh @@ -3,10 +3,10 @@ #include "Python.h" -#include "ExternalEffects.hh" +#include #include "BaseLaserFieldSource.hh" #include "Stark.hh" -#include "ParticleAttributes.hh" +#include using namespace TrackerRK4; diff --git a/DiagonalDensityMatrix/wrap_DM_noLaserField.cc b/DiagonalDensityMatrix/wrap_DM_noLaserField.cc index 947ccc0..49672da 100644 --- a/DiagonalDensityMatrix/wrap_DM_noLaserField.cc +++ b/DiagonalDensityMatrix/wrap_DM_noLaserField.cc @@ -1,7 +1,6 @@ -#include "orbit_mpi.hh" -#include "pyORBIT_Object.hh" +#include +#include -#include "wrap_utils.hh" #include "wrap_DM_noLaserField.hh" #include @@ -40,28 +39,28 @@ extern "C" { //this is implementation of the __init__ method static int DM_noLaserField_init(pyORBIT_Object *self, PyObject *args, PyObject *kwds){ PyObject* pyStarkEffect=NULL; - + int nVars = PyTuple_Size(args); int method = 0; - + if(!PyArg_ParseTuple( args,"Oi:",&pyStarkEffect, &method)) { error("DM_noLaserField(StarkEffect) - param. needed"); - } - + } + else { - + Stark* Starkef = (Stark*) ((pyORBIT_Object*) pyStarkEffect)->cpp_obj; self->cpp_obj = new DM_noLaserField(Starkef, method); ((DM_noLaserField*) self->cpp_obj)->setPyWrapper((PyObject*) self); } - + return 0; } - - - - // name([name]) - sets or returns the name of the External Effeects class + + + + // name([name]) - sets or returns the name of the External Effeects class static PyObject* DM_noLaserField_name(PyObject *self, PyObject *args){ DM_noLaserField* cpp_DM_noLaserField = (DM_noLaserField*) ((pyORBIT_Object*) self)->cpp_obj; const char* name = NULL; @@ -73,15 +72,15 @@ extern "C" { cpp_DM_noLaserField->setName(name_str); } return Py_BuildValue("s",cpp_DM_noLaserField->getName().c_str()); - } - + } + //----------------------------------------------------- //destructor for python PyExternalEffects class (__del__ method). //----------------------------------------------------- static void DM_noLaserField_del(pyORBIT_Object* self){ //std::cerr<<"The DM_noLaserField __del__ has been called!"<cpp_obj); - self->ob_type->tp_free((PyObject*)self); + Py_TYPE(self)->tp_free((PyObject*)self); } // defenition of the methods of the python PyExternalEffects wrapper class @@ -99,8 +98,7 @@ extern "C" { //new python PyExternalEffects wrapper type definition static PyTypeObject pyORBIT_DM_noLaserField_Type = { - PyObject_HEAD_INIT(NULL) - 0, /*ob_size*/ + PyVarObject_HEAD_INIT(NULL, 0) "DM_noLaserField", /*tp_name*/ sizeof(pyORBIT_Object), /*tp_basicsize*/ 0, /*tp_itemsize*/ @@ -138,10 +136,10 @@ extern "C" { (initproc) DM_noLaserField_init, /* tp_init */ 0, /* tp_alloc */ DM_noLaserField_new, /* tp_new */ - }; + }; + - //-------------------------------------------------- //Initialization function of the pyPyExternalEffects class //It will be called from Bunch wrapper initialization @@ -150,7 +148,7 @@ extern "C" { if (PyType_Ready(&pyORBIT_DM_noLaserField_Type) < 0) return; Py_INCREF(&pyORBIT_DM_noLaserField_Type); PyModule_AddObject(module, "DM_noLaserField", (PyObject *)&pyORBIT_DM_noLaserField_Type); - + } #ifdef __cplusplus diff --git a/FringeField/FringeField.cc b/FringeField/FringeField.cc index 9ffcdf9..41bcd5b 100644 --- a/FringeField/FringeField.cc +++ b/FringeField/FringeField.cc @@ -7,8 +7,8 @@ // 04/21/2003 // // DESCRIPTION -// The base class for Python implementation of a field source. -// It should be sub-classed on Python level and implements +// The base class for Python implementation of a field source. +// It should be sub-classed on Python level and implements // getElectricField(x,y,z,t) and getMagneticField (x,y,z,t) methods. // The results of these methods will be available from the c++ level. // This is an example of embedding Python in C++ Orbit level. @@ -17,8 +17,8 @@ /////////////////////////////////////////////////////////////////////////// #include "FringeField.hh" -#include "orbit_mpi.hh" -#include "BufferStore.hh" +#include +#include #include #include @@ -28,7 +28,7 @@ using namespace OrbitUtils; FringeField::FringeField(double d,double F, int derivative) { - + gap = d; Bx = F; der = derivative; @@ -39,7 +39,7 @@ FringeField::FringeField(double d,double F, int derivative) FringeField::~FringeField() -{ +{ } @@ -50,27 +50,27 @@ FringeField::~FringeField() -void FringeField::getElectricMagneticField(double x, double y, double z, double t, +void FringeField::getElectricMagneticField(double x, double y, double z, double t, double& E_x, double& E_y, double& E_z, double& H_x, double& H_y, double& H_z) { double coeff = -der*z/gap; - + double coeff_2 = coeff; double coeff_3 = coeff*coeff_2; double coeff_4 = coeff*coeff_3; double coeff_5 = coeff*coeff_4; double coeff_6 = coeff*coeff_5; - - + + E_x=0; E_y=0; E_z=0; - + H_y=0; H_z=0; - + H_x=Bx/(1+exp(0.478959+1.911289*coeff_2-1.185953*coeff_3+1.630554*coeff_4-1.082657*coeff_5+0.318111*coeff_6)); diff --git a/FringeField/FringeField.hh b/FringeField/FringeField.hh index 10187d5..48adb51 100644 --- a/FringeField/FringeField.hh +++ b/FringeField/FringeField.hh @@ -23,7 +23,7 @@ #define FRINGEFIELD_HH_ #include "Python.h" -#include "BaseFieldSource.hh" +#include #include "FieldOrientation.hh" #include diff --git a/FringeField/wrap_fringe_field.cc b/FringeField/wrap_fringe_field.cc index b1b897e..3bd2f2c 100644 --- a/FringeField/wrap_fringe_field.cc +++ b/FringeField/wrap_fringe_field.cc @@ -1,7 +1,6 @@ -#include "orbit_mpi.hh" -#include "pyORBIT_Object.hh" +#include +#include -#include "wrap_utils.hh" #include "wrap_fringe_field.hh" #include @@ -10,7 +9,6 @@ #include "FringeField.hh" using namespace OrbitUtils; -using namespace wrap_orbit_utils; namespace wrap_fringe_field{ @@ -37,9 +35,9 @@ extern "C" { //initializator for python CppBaseFieldSource class //this is implementation of the __init__ method static int FringeField_init(pyORBIT_Object *self, PyObject *args, PyObject *kwds){ - - + + double d; double F; @@ -47,28 +45,28 @@ extern "C" { if(!PyArg_ParseTuple( args,"ddi:",&d, &F, &der)){ error("FringeField(gap , By_max, derivative of module |By(z)| (plus one or minus one)) - params. are needed"); - } + } else { self->cpp_obj = new FringeField(d, F, der); ((FringeField*) self->cpp_obj)->setPyWrapper((PyObject*) self); } - - + + return 0; } - - - - - + + + + + static PyObject* FringeField_getField(PyObject *self, PyObject *args){ FringeField* cpp_fields = (FringeField*)((pyORBIT_Object*) self)->cpp_obj; - - + + int nVars = PyTuple_Size(args); - + double x; double y; double z; @@ -79,7 +77,7 @@ extern "C" { double Bx; double By; double Bz; - + //NO NEW OBJECT CREATED BY PyArg_ParseTuple! NO NEED OF Py_DECREF() if(!PyArg_ParseTuple( args,"dddd:",&x,&y,&z,&t)) @@ -89,8 +87,8 @@ extern "C" { return Py_BuildValue("d",Bx); } - - + + //----------------------------------------------------- @@ -99,7 +97,7 @@ extern "C" { static void FringeField_del(pyORBIT_Object* self){ //std::cerr<<"The FringeField __del__ has been called!"<cpp_obj); - self->ob_type->tp_free((PyObject*)self); + Py_TYPE(self)->tp_free((PyObject*)self); } // defenition of the methods of the python PyBaseFieldSource wrapper class @@ -125,8 +123,7 @@ extern "C" { //new python PyBaseFieldSource wrapper type definition static PyTypeObject pyORBIT_FringeField_Type = { - PyObject_HEAD_INIT(NULL) - 0, /*ob_size*/ + PyVarObject_HEAD_INIT(NULL, 0) "FringeField", /*tp_name*/ sizeof(pyORBIT_Object), /*tp_basicsize*/ 0, /*tp_itemsize*/ @@ -164,7 +161,7 @@ extern "C" { (initproc) FringeField_init, /* tp_init */ 0, /* tp_alloc */ FringeField_new, /* tp_new */ - }; + }; //-------------------------------------------------- //Initialization function of the pyPyBaseFieldSource class diff --git a/FroissartStoraLF/FroissartStoraLF.cc b/FroissartStoraLF/FroissartStoraLF.cc index 29b25ed..fe4d3f7 100644 --- a/FroissartStoraLF/FroissartStoraLF.cc +++ b/FroissartStoraLF/FroissartStoraLF.cc @@ -7,8 +7,8 @@ // 04/21/2003 // // DESCRIPTION -// The base class for Python implementation of a field source. -// It should be sub-classed on Python level and implements +// The base class for Python implementation of a field source. +// It should be sub-classed on Python level and implements // getElectricField(x,y,z,t) and getMagneticField (x,y,z,t) methods. // The results of these methods will be available from the c++ level. // This is an example of embedding Python in C++ Orbit level. @@ -16,7 +16,7 @@ // /////////////////////////////////////////////////////////////////////////// #include "FroissartStoraLF.hh" -#include "orbit_mpi.hh" +#include #include #define J tcomplex(0.,1.) @@ -25,17 +25,17 @@ using namespace OrbitUtils; FroissartStoraLF::FroissartStoraLF(double Om,double G,double El) { - + Omega=Om; Gamma=G; amp_ELas = El; - + //default polarization of laser field nEx=1; nEy=0; nEz=0; - + } @@ -53,14 +53,14 @@ FroissartStoraLF::~FroissartStoraLF() void FroissartStoraLF::setLaserFieldPolarization(double n_Ex, double n_Ey, double n_Ez) { - + double a=sqrt(n_Ex*n_Ex+n_Ey*n_Ey+n_Ez*n_Ez); - - nEx=n_Ex/a; - nEy=n_Ey/a; + + nEx=n_Ex/a; + nEy=n_Ey/a; nEz=n_Ez/a; //Vector of polarization of electric field - + } @@ -71,14 +71,14 @@ void FroissartStoraLF::setLaserFieldPolarization(double n_Ex, double n_Ey, doubl tcomplex FroissartStoraLF::getNonOrientedU(double t){ - + double ta=2.418884326505e-17; //atomic unit of time double Ea=5.14220642e11; //Atomic unit of electric field t/=ta; //convertion of natural time to atomic unit time - + return amp_ELas*Ea*exp(J*(Omega*t+Gamma*t*t/2)); - + } @@ -89,19 +89,19 @@ tcomplex FroissartStoraLF::getNonOrientedU(double t){ -tcomplex FroissartStoraLF::getLaserEMField(double x, double y, double z, double t, +tcomplex FroissartStoraLF::getLaserEMField(double x, double y, double z, double t, tcomplex& E_x, tcomplex& E_y, tcomplex& E_z, tcomplex& H_x, tcomplex& H_y, tcomplex& H_z){ - - + + tcomplex E=getNonOrientedU(t); E_x=E*nEx; E_y=E*nEy; E_z=E*nEz; H_x=0; H_y=0; H_z=0; - + return E/abs(E); - + } @@ -113,9 +113,9 @@ double FroissartStoraLF::getFrequencyOmega(double m, double x, double y, double double ta=2.418884326505e-17; //atomic unit of time t/=ta; //convertion of natural time to atomic unit time - + return (Omega+Gamma*t)/ta; - + } diff --git a/FroissartStoraLF/wrap_froissart_stora_lf.cc b/FroissartStoraLF/wrap_froissart_stora_lf.cc index ce315ad..371d6d1 100644 --- a/FroissartStoraLF/wrap_froissart_stora_lf.cc +++ b/FroissartStoraLF/wrap_froissart_stora_lf.cc @@ -1,7 +1,6 @@ -#include "orbit_mpi.hh" -#include "pyORBIT_Object.hh" +#include +#include -#include "wrap_utils.hh" #include "wrap_froissart_stora_lf.hh" #include @@ -12,7 +11,6 @@ #include "FroissartStoraLF.hh" using namespace OrbitUtils; -using namespace wrap_orbit_utils; namespace wrap_froissart_stora_lf{ @@ -39,31 +37,31 @@ extern "C" { //initializator for python CppBaseFieldSource class //this is implementation of the __init__ method static int FroissartStoraLF_init(pyORBIT_Object *self, PyObject *args, PyObject *kwds){ - + double Omega; double Gamma; double Elas; - + if(!PyArg_ParseTuple( args,"ddd:",&Omega,&Gamma,&Elas)){ error("FroissartStoraLF(Omega,Gamma,ampl_Elas) - params. are needed"); - } + } else { - + self->cpp_obj = new FroissartStoraLF(Omega,Gamma,Elas); ((FroissartStoraLF*) self->cpp_obj)->setPyWrapper((PyObject*) self); } return 0; } - - - - - - + + + + + + static PyObject* FroissartStoraLF_setLaserFieldPolarization(PyObject *self, PyObject *args){ FroissartStoraLF* LaserField = (FroissartStoraLF*)((pyORBIT_Object*) self)->cpp_obj; - + double nEx; double nEy; @@ -72,19 +70,19 @@ extern "C" { if(!PyArg_ParseTuple( args,"ddd:", &nEx, &nEy, &nEz)) {error("LaserExternalEfects - setLaserFieldPolarization(nEx, nEy, nEz) - params. afe needed");} - - else + + else LaserField->setLaserFieldPolarization(nEx, nEy, nEz); - - + + Py_INCREF(Py_None); - return Py_None; - } - - - - - + return Py_None; + } + + + + + //----------------------------------------------------- //destructor for python PyBaseFieldSource class (__del__ method). @@ -92,7 +90,7 @@ extern "C" { static void FroissartStoraLF_del(pyORBIT_Object* self){ //std::cerr<<"The CppBaseFieldSource __del__ has been called!"<cpp_obj); - self->ob_type->tp_free((PyObject*)self); + Py_TYPE(self)->tp_free((PyObject*)self); } // defenition of the methods of the python PyBaseFieldSource wrapper class @@ -118,8 +116,7 @@ extern "C" { //new python PyBaseFieldSource wrapper type definition static PyTypeObject pyORBIT_FroissartStoraLF_Type = { - PyObject_HEAD_INIT(NULL) - 0, /*ob_size*/ + PyVarObject_HEAD_INIT(NULL, 0) "FroissartStoraLF", /*tp_name*/ sizeof(pyORBIT_Object), /*tp_basicsize*/ 0, /*tp_itemsize*/ @@ -157,7 +154,7 @@ extern "C" { (initproc) FroissartStoraLF_init, /* tp_init */ 0, /* tp_alloc */ FroissartStoraLF_new, /* tp_new */ - }; + }; //-------------------------------------------------- //Initialization function of the pyPyBaseFieldSource class diff --git a/HermiteGaussianLFmode/HermiteGaussianLFmode.cc b/HermiteGaussianLFmode/HermiteGaussianLFmode.cc index c528083..abfe4d3 100644 --- a/HermiteGaussianLFmode/HermiteGaussianLFmode.cc +++ b/HermiteGaussianLFmode/HermiteGaussianLFmode.cc @@ -7,8 +7,8 @@ // 04/21/2003 // // DESCRIPTION -// The base class for Python implementation of a field source. -// It should be sub-classed on Python level and implements +// The base class for Python implementation of a field source. +// It should be sub-classed on Python level and implements // getElectricField(x,y,z,t) and getMagneticField (x,y,z,t) methods. // The results of these methods will be available from the c++ level. // This is an example of embedding Python in C++ Orbit level. @@ -16,10 +16,10 @@ // /////////////////////////////////////////////////////////////////////////// #include "HermiteGaussianLFmode.hh" -#include "MathPolynomial.hh" +#include -#include "orbit_mpi.hh" -#include "OrbitConst.hh" +#include +#include #include @@ -32,7 +32,7 @@ HermiteGaussianLFmode::HermiteGaussianLFmode(tcomplex C,int n,int m,double w_x,d Cnm=C; Unm=Cnm*sqrt(2*OrbitConst::c*OrbitConst::permeability)*sqrt(pow(w_x,2*n+1)*pow(w_y,2*m+1)*MathPolynomial::Factorial(n)*MathPolynomial::Factorial(m)*pow(2,n+m+1)/(MathPolynomial::Factorial(2*n)*MathPolynomial::Factorial(2*m)*OrbitConst::PI)); - orient = new FieldOrientation(); + orient = new FieldOrientation(); Laser_lambda=la; env_peak = env__peak; env_sigma = env__sigma; @@ -43,7 +43,7 @@ HermiteGaussianLFmode::HermiteGaussianLFmode(tcomplex C,int n,int m,double w_x,d n_moda=n; m_moda=m; k=2*OrbitConst::PI/Laser_lambda; - + la_pi_wx = Laser_lambda/(OrbitConst::PI*wx); la_pi_wy = Laser_lambda/(OrbitConst::PI*wy); @@ -52,13 +52,13 @@ HermiteGaussianLFmode::HermiteGaussianLFmode(tcomplex C,int n,int m,double w_x,d nHx=0;nHy=1;nHz=0; - - - + + + } HermiteGaussianLFmode::~HermiteGaussianLFmode() -{ +{ } @@ -72,26 +72,26 @@ HermiteGaussianLFmode::~HermiteGaussianLFmode() void HermiteGaussianLFmode::setLaserFieldOrientation(double x0, double y0, double z0, double kx, double ky, double kz, double mx, double my, double mz, - double n_Ex, double n_Ey, double n_Ez) + double n_Ex, double n_Ey, double n_Ez) { double a=sqrt(n_Ex*n_Ex+n_Ey*n_Ey+n_Ez*n_Ez); orient->setCoefficients(x0, y0, z0,kx, ky, kz, mx, my, mz); - - nEx=n_Ex/a; - nEy=n_Ey/a; + + nEx=n_Ex/a; + nEy=n_Ey/a; nEz=n_Ez/a; //Vector of polarization of electric field - + double n_Hx=nEz*ky-nEy*kz; double n_Hy=nEx*kz-nEz*kx; double n_Hz=nEy*kx-nEx*ky; - + double b=sqrt(n_Hx*n_Hx+n_Hy*n_Hy+n_Hz*n_Hz); - + nHx=n_Hx/b; nHy=n_Hy/b; nHz=n_Hz/b; - + } @@ -99,17 +99,17 @@ double b=sqrt(n_Hx*n_Hx+n_Hy*n_Hy+n_Hz*n_Hz); void HermiteGaussianLFmode::setLocalParameters(double rx, double ry, double ax, double ay) { - - + + fx=-4*k*k*ax*rx*rx*rx/(1+4*ax*ax*rx*rx*k*k); fy=-4*k*k*ay*ry*ry*ry/(1+4*ay*ay*ry*ry*k*k); - + wx=2*rx/sqrt(1+4*ax*ax*rx*rx*k*k); wy=2*ry/sqrt(1+4*ay*ay*ry*ry*k*k); - + int n=n_moda; int m=m_moda; - + Unm=Cnm*sqrt(2*OrbitConst::c*OrbitConst::permeability)*sqrt(pow(wx,2*n+1)*pow(wy,2*m+1)*MathPolynomial::Factorial(n)*MathPolynomial::Factorial(m)*pow(2,n+m+1)/(MathPolynomial::Factorial(2*n)*MathPolynomial::Factorial(2*m)*OrbitConst::PI)); @@ -119,15 +119,15 @@ void HermiteGaussianLFmode::setLocalParameters(double rx, double ry, double ax, double HermiteGaussianLFmode::HermiteGaussianOmega(double m,double x, double y, double z, double px, double py, double pz, double t){ - - - + + + double a=OrbitConst::c/sqrt(m*m+px*px+py*py+pz*pz); - + double vx=px*a; double vy=py*a; double vz=pz*a; - + double ax=k*k*pow(wx,4)+4*pow(fx-z,2); @@ -151,18 +151,18 @@ double by = 4*y*(fy-z); -tcomplex HermiteGaussianLFmode::getLaserEMField(double x, double y, double z, double t, +tcomplex HermiteGaussianLFmode::getLaserEMField(double x, double y, double z, double t, tcomplex& E_x, tcomplex& E_y, tcomplex& E_z, tcomplex& H_x, tcomplex& H_y, tcomplex& H_z){ - + tcomplex exp_phasa; orient->OrientCoordinates(x,y,z); - + tcomplex E=Unm*getNonOrientedU(n_moda,m_moda,x,y,z,t); exp_phasa = E/abs(E); - double absp = abs(exp_phasa); + double absp = abs(exp_phasa); if ((absp != absp)||(fabs(absp-1.)>1.0e-3)) exp_phasa = tcomplex(0.,0.); @@ -177,22 +177,22 @@ tcomplex H=E/OrbitConst::c; E_x=E*nEx; E_y=E*nEy; E_z=E*nEz; H_x=H*nHx; H_y=H*nHy; H_z=H*nHz; - + return exp_phasa; - + } bool HermiteGaussianLFmode::region(double x, double y, double z){ - + double ax = (z-fx)*la_pi_wx; double ay = (z-fy)*la_pi_wy; - - + + orient->OrientCoordinates(x,y,z); - + return (x*x/(wx*wx+ax*ax) + y*y/(wy*wy+ay*ay) < 25.); - + } @@ -200,10 +200,10 @@ bool HermiteGaussianLFmode::region(double x, double y, double z){ double HermiteGaussianLFmode::getFrequencyOmega(double m, double x, double y, double z, double px, double py, double pz, double t){ - + orient->OrientCoordinates(x,y,z); orient->OrientVector2(px,py,pz); - + return HermiteGaussianOmega(m,x,y,z,px,py,pz,t); @@ -222,36 +222,36 @@ double HermiteGaussianLFmode::getFrequencyOmega(double m, double x, double y, do tcomplex HermiteGaussianLFmode::getNonOrientedU(int n, int m, double x, double y, double z, double t){ tcomplex U; - + if (n==0&&m==0) { - + tcomplex a=tcomplex(0.,k*(z-t*OrbitConst::c)); tcomplex funx=pow(tcomplex(wx*wx,-(z-fx)*2/k),-1); tcomplex funy=pow(tcomplex(wy*wy,-(z-fy)*2/k),-1); - - U = sqrt(funx*funy)*exp(-x*x*funx-y*y*funy-a); + + U = sqrt(funx*funy)*exp(-x*x*funx-y*y*funy-a); if (abs(U) != abs(U)) U = tcomplex(0.,0.); - - return U; - + + return U; + } else { - - - + + + tcomplex funx=sqrt(tcomplex(wx*wx-2.0*J*(z-fx)/k)); tcomplex funy=sqrt(tcomplex(wy*wy-2.0*J*(z-fy)/k)); tcomplex xf=x/funx; tcomplex yf=y/funy; - + U = pow(funx,-n-1)*pow(funy,-m-1)*exp(-xf*xf-yf*yf-J*k*(z-t*OrbitConst::c))*MathPolynomial::ComplexHermite(n,xf)*MathPolynomial::ComplexHermite(m,yf); if (abs(U) != abs(U)) U = tcomplex(0.,0.); - + return U; - + } - - + + } diff --git a/HermiteGaussianLFmode/wrap_hermite_gaussian_lf_mode.cc b/HermiteGaussianLFmode/wrap_hermite_gaussian_lf_mode.cc index a2be284..cf5032d 100644 --- a/HermiteGaussianLFmode/wrap_hermite_gaussian_lf_mode.cc +++ b/HermiteGaussianLFmode/wrap_hermite_gaussian_lf_mode.cc @@ -1,7 +1,6 @@ -#include "orbit_mpi.hh" -#include "pyORBIT_Object.hh" +#include +#include -#include "wrap_utils.hh" #include "wrap_hermite_gaussian_lf_mode.hh" #include @@ -12,7 +11,6 @@ #include "HermiteGaussianLFmode.hh" using namespace OrbitUtils; -using namespace wrap_orbit_utils; namespace wrap_hermite_gaussian_lf_mode{ @@ -39,7 +37,7 @@ extern "C" { //initializator for python CppBaseFieldSource class //this is implementation of the __init__ method static int HermiteGaussianLFmode_init(pyORBIT_Object *self, PyObject *args, PyObject *kwds){ - + double Cnm=0; int n=0; int m=0; @@ -61,14 +59,14 @@ extern "C" { return 0; } - - - - - - + + + + + + static PyObject* HermiteGaussianLFmode_setLaserFieldOrientation(PyObject *self, PyObject *args){ - + HermiteGaussianLFmode* LaserField = (HermiteGaussianLFmode*)((pyORBIT_Object*) self)->cpp_obj; double x0; double y0; @@ -86,25 +84,25 @@ extern "C" { if(!PyArg_ParseTuple( args,"dddddddddddd:",&x0, &y0, &z0, &kx, &ky, &kz, &mx, &my, &mz, &nEx, &nEy, &nEz)) {error("LaserExternalEfects - setLaserHalfAngle(comments is not written yet) - param. is needed");flag=false;} - + if(fabs(kx*mx+ky*my+kz*mz)/sqrt(kx*kx+ky*ky+kz*kz)/sqrt(mx*mx+my*my+mz*mz)>1.e-15) {error("Please be shure that kx*mx+ky*my+kz*mz==0");flag=false;} if(fabs(kx*nEx+ky*nEy+kz*nEz)/sqrt(kx*kx+ky*ky+kz*kz)/sqrt(nEx*nEx+nEy*nEy+nEz*nEz)>1.e-15) {error("Please be shure that kx*Ex+ky*Ey+kz*Ez==0");flag=false;} - - if(flag) + + if(flag) LaserField->setLaserFieldOrientation(x0,y0, z0, kx, ky, kz, mx, my, mz, nEx, nEy, nEz); - - + + Py_INCREF(Py_None); - return Py_None; - } - - - + return Py_None; + } + + + static PyObject* HermiteGaussianLFmode_setLocalParameters(PyObject *self, PyObject *args){ - + HermiteGaussianLFmode* LaserField = (HermiteGaussianLFmode*)((pyORBIT_Object*) self)->cpp_obj; double rx; @@ -117,16 +115,16 @@ extern "C" { error("LaserExternalEfects - setLaserHalfAngle(comments is not written yet) - param. is needed"); else LaserField->setLocalParameters(rx,ry,ax,ay); - - - + + + Py_INCREF(Py_None); - return Py_None; - } - - - - + return Py_None; + } + + + + /* void LasStripExternalEffects::setLaserHalfAngle(double a); double LasStripExternalEffects::getLaserHalfAngle(); @@ -135,59 +133,59 @@ extern "C" { void LasStripExternalEffects::setLaser_lambda(double a); double LasStripExternalEffects::getLaser_lambda(); - + static PyObject* HermiteGaussianLFmode_setLaserHalfAngle(PyObject *self, PyObject *args){ HermiteGaussianLFmode* LaserField = (HermiteGaussianLFmode*)((pyORBIT_Object*) self)->cpp_obj; - double half_Angle; + double half_Angle; if(!PyArg_ParseTuple( args,"d:",&half_Angle)){ error("LaserExternalEfects - setLaserHalfAngle(half_Angle) - param. is needed"); } else LaserField->setLaserHalfAngle(half_Angle); - - + + Py_INCREF(Py_None); - return Py_None; - } + return Py_None; + } static PyObject* HermiteGaussianLFmode_setLaserPower(PyObject *self, PyObject *args){ HermiteGaussianLFmode* LaserField = (HermiteGaussianLFmode*)((pyORBIT_Object*) self)->cpp_obj; - double LaserPower; + double LaserPower; if(!PyArg_ParseTuple( args,"d:",&LaserPower)){ error("LaserExternalEfects - setLaserPower(LaserPower) - param. is needed"); } else LaserField->setLaserPower(LaserPower); - - + + Py_INCREF(Py_None); - return Py_None; + return Py_None; } static PyObject* HermiteGaussianLFmode_setLaser_lambda(PyObject *self, PyObject *args){ HermiteGaussianLFmode* LaserField = (HermiteGaussianLFmode*)((pyORBIT_Object*) self)->cpp_obj; - double Laser_lambda; + double Laser_lambda; if(!PyArg_ParseTuple( args,"d:",&Laser_lambda)){ error("LaserExternalEfects - setLaser_lambda(Laser_lambda) - param. is needed"); } else LaserField->setLaserPower(Laser_lambda); - - + + Py_INCREF(Py_None); - return Py_None; + return Py_None; } - - - + + + static PyObject* HermiteGaussianLFmode_setLaserFieldOrientation(PyObject *self, PyObject *args){ HermiteGaussianLFmode* LaserField = (HermiteGaussianLFmode*)((pyORBIT_Object*) self)->cpp_obj; - double Laser_lambda; + double Laser_lambda; if(!PyArg_ParseTuple( args,"d:",&Laser_lambda)){ error("LaserExternalEfects - setLaser_lambda(Laser_lambda) - param. is needed"); } else LaserField->setLaserPower(Laser_lambda); - - + + Py_INCREF(Py_None); - return Py_None; + return Py_None; } @@ -195,28 +193,28 @@ extern "C" { static PyObject* HermiteGaussianLFmode_getLaserHalfAngle(PyObject *self, PyObject *args){ HermiteGaussianLFmode* LaserField = (HermiteGaussianLFmode*)((pyORBIT_Object*) self)->cpp_obj; - - return Py_BuildValue("d",LaserField->getLaserHalfAngle()); - } - + + return Py_BuildValue("d",LaserField->getLaserHalfAngle()); + } + static PyObject* HermiteGaussianLFmode_getLaserPower(PyObject *self, PyObject *args){ HermiteGaussianLFmode* LaserField = (HermiteGaussianLFmode*)((pyORBIT_Object*) self)->cpp_obj; - - return Py_BuildValue("d",LaserField->getLaserPower()); - } + + return Py_BuildValue("d",LaserField->getLaserPower()); + } static PyObject* HermiteGaussianLFmode_getLaser_lambda(PyObject *self, PyObject *args){ HermiteGaussianLFmode* LaserField = (HermiteGaussianLFmode*)((pyORBIT_Object*) self)->cpp_obj; - - return Py_BuildValue("d",LaserField->getLaser_omega()); - } + + return Py_BuildValue("d",LaserField->getLaser_omega()); + } */ - - - + + + //----------------------------------------------------- //destructor for python PyBaseFieldSource class (__del__ method). @@ -224,7 +222,7 @@ extern "C" { static void HermiteGaussianLFmode_del(pyORBIT_Object* self){ //std::cerr<<"The CppBaseFieldSource __del__ has been called!"<cpp_obj); - self->ob_type->tp_free((PyObject*)self); + Py_TYPE(self)->tp_free((PyObject*)self); } // defenition of the methods of the python PyBaseFieldSource wrapper class @@ -250,8 +248,7 @@ extern "C" { //new python PyBaseFieldSource wrapper type definition static PyTypeObject pyORBIT_HermiteGaussianLFmode_Type = { - PyObject_HEAD_INIT(NULL) - 0, /*ob_size*/ + PyVarObject_HEAD_INIT(NULL, 0) "HermiteGaussianLFmode", /*tp_name*/ sizeof(pyORBIT_Object), /*tp_basicsize*/ 0, /*tp_itemsize*/ @@ -289,7 +286,7 @@ extern "C" { (initproc) HermiteGaussianLFmode_init, /* tp_init */ 0, /* tp_alloc */ HermiteGaussianLFmode_new, /* tp_new */ - }; + }; //-------------------------------------------------- //Initialization function of the pyPyBaseFieldSource class diff --git a/HminusStripping/HminusStripping.cc b/HminusStripping/HminusStripping.cc index 821f805..8e56626 100644 --- a/HminusStripping/HminusStripping.cc +++ b/HminusStripping/HminusStripping.cc @@ -7,26 +7,26 @@ // 06/27/2008 // // DESCRIPTION -// The base class for C++ implementation of a external effects -// during the transport of particles through the external field. +// The base class for C++ implementation of a external effects +// during the transport of particles through the external field. // It should be sub-classed on Python level and implement // setupEffects(Bunch* bunch) // finalizeEffects(Bunch* bunch) -// applyEffects(Bunch* bunch, int index, -// double* y_in_vct, double* y_out_vct, -// double t, double t_step, +// applyEffects(Bunch* bunch, int index, +// double* y_in_vct, double* y_out_vct, +// double t, double t_step, // OrbitUtils::BaseFieldSource* fieldSource) // methods. // The results of these methods will be available from the c++ level. // This is an example of embedding Python in C++ Orbit level. // -// -// -// -// +// +// +// +// // /////////////////////////////////////////////////////////////////////////// -#include "orbit_mpi.hh" +#include #include #include #include @@ -36,9 +36,9 @@ #include "HminusStripping.hh" -#include "RungeKuttaTracker.hh" -#include "OrbitConst.hh" -#include "LorentzTransformationEM.hh" +#include +#include +//#include "LorentzTransformationEM.hh" @@ -63,8 +63,8 @@ HminusStripping::HminusStripping(int method) setName("unnamed"); index = method; - - + + //allocating memory for koefficients of 4-th order Runge-Kutta method and other koeeficients of the master equation @@ -76,7 +76,7 @@ k_RungeKutt=new double[5]; HminusStripping::~HminusStripping() { - + delete [] k_RungeKutt; } @@ -85,10 +85,10 @@ HminusStripping::~HminusStripping() -void HminusStripping::setupEffects(Bunch* bunch){ - +void HminusStripping::setupEffects(Bunch* bunch){ + + - prob = new double[bunch->getSize()]; bunch->setCharge(-1); @@ -98,7 +98,7 @@ void HminusStripping::setupEffects(Bunch* bunch){ //std::cout<hasParticleAttributes("Populations")!=0) if(bunch->getParticleAttributes("Populations")->getAttSize() != 1) @@ -113,20 +113,20 @@ void HminusStripping::setupEffects(Bunch* bunch){ for (int i=0; igetSize();i++) bunch->getParticleAttributes("Populations")->attValue(i,0) = 0.; } - + PopAttr = bunch->getParticleAttributes("Populations"); - + } - - + + void HminusStripping::finalizeEffects(Bunch* bunch){ - + delete [] prob; - + } @@ -136,25 +136,25 @@ void HminusStripping::finalizeEffects(Bunch* bunch){ -void HminusStripping::applyEffectsForEach(Bunch* bunch, int i, - double* y_in_vct, double* y_out_vct, +void HminusStripping::applyEffectsForEach(Bunch* bunch, int i, + double* y_in_vct, double* y_out_vct, double t, double t_step, BaseFieldSource* fieldSource, RungeKuttaTracker* tracker) { - + if(bunch->flag(i) == 0) { - + double** partCoordArr = bunch->coordArr(); - + y_in_vct[0] = partCoordArr[i][0]; y_in_vct[1] = partCoordArr[i][2]; y_in_vct[2] = partCoordArr[i][4]; y_in_vct[3] = partCoordArr[i][1]; y_in_vct[4] = partCoordArr[i][3]; - y_in_vct[5] = partCoordArr[i][5]; - - double mass = bunch->getMass(); + y_in_vct[5] = partCoordArr[i][5]; + + double mass = bunch->getMass(); double vpt = mass*mass + y_in_vct[3]*y_in_vct[3] + y_in_vct[4]*y_in_vct[4] + y_in_vct[5]*y_in_vct[5]; - + vpt = 299792458*t_step/sqrt(vpt); y_out_vct[0] = y_in_vct[0] + vpt*y_in_vct[3]; y_out_vct[1] = y_in_vct[1] + vpt*y_in_vct[4]; @@ -162,41 +162,41 @@ void HminusStripping::applyEffectsForEach(Bunch* bunch, int i, y_out_vct[3] = y_in_vct[3]; y_out_vct[4] = y_in_vct[4]; y_out_vct[5] = y_in_vct[5]; - + partCoordArr[i][0] = y_out_vct[0]; partCoordArr[i][2] = y_out_vct[1]; partCoordArr[i][4] = y_out_vct[2]; partCoordArr[i][1] = y_out_vct[3]; partCoordArr[i][3] = y_out_vct[4]; - partCoordArr[i][5] = y_out_vct[5]; - + partCoordArr[i][5] = y_out_vct[5]; + } - + x0 = y_in_vct[0];y0 = y_in_vct[1];z0 = y_in_vct[2]; px0 = y_in_vct[3];py0 = y_in_vct[4];pz0 = y_in_vct[5]; x = y_out_vct[0];y = y_out_vct[1];z = y_out_vct[2]; px = y_out_vct[3];py = y_out_vct[4];pz = y_out_vct[5]; - - - // This function gives parameters Ez_stat Ex_las[1...3] Ey_las[1...3] Ez_las[1...3] - //in natural unts (Volt per meter) in the frame of particle + + + // This function gives parameters Ez_stat Ex_las[1...3] Ey_las[1...3] Ez_las[1...3] + //in natural unts (Volt per meter) in the frame of particle GetParticleFrameFields(i, t,t_step, bunch,fieldSource); - - // This function gives parameters Ez_stat Ex_las[1...3] Ey_las[1...3] Ez_las[1...3] t_part omega_part part_t_step - //in atomic units in frame of particle - GetParticleFrameParameters(i,t,t_step,bunch); + + // This function gives parameters Ez_stat Ex_las[1...3] Ey_las[1...3] Ez_las[1...3] t_part omega_part part_t_step + //in atomic units in frame of particle + GetParticleFrameParameters(i,t,t_step,bunch); // This function provides step solution for density matrix using rk4 method AmplSolver4step(t_step, i,bunch); - - + + // cout<x(0)<<"\t"<y(0)<<"\t"<z(0)<<"\n"; - + } @@ -215,24 +215,24 @@ void HminusStripping::applyEffectsForEach(Bunch* bunch, int i, void HminusStripping::GetParticleFrameFields(int i,double t,double t_step, Bunch* bunch, BaseFieldSource* fieldSource) { - + double Ez; - - fieldSource->getElectricMagneticField(x0,y0,z0,t,Ex_stat,Ey_stat,Ez_stat,Bx_stat,By_stat,Bz_stat); - LorentzTransformationEM::transform(bunch->getMass(),px0,py0,pz0,Ex_stat,Ey_stat,Ez_stat,Bx_stat,By_stat,Bz_stat); - - - + fieldSource->getElectricMagneticField(x0,y0,z0,t,Ex_stat,Ey_stat,Ez_stat,Bx_stat,By_stat,Bz_stat); + //LorentzTransformationEM::transform(bunch->getMass(),px0,py0,pz0,Ex_stat,Ey_stat,Ez_stat,Bx_stat,By_stat,Bz_stat); + + + + Ez=sqrt(Ex_stat*Ex_stat+Ey_stat*Ey_stat+Ez_stat*Ez_stat); - - Ex_stat=0; - Ey_stat=0; + + Ex_stat=0; + Ey_stat=0; Ez_stat=Ez; - - - + + + } @@ -243,8 +243,8 @@ Ez=sqrt(Ex_stat*Ex_stat+Ey_stat*Ey_stat+Ez_stat*Ez_stat); void HminusStripping::GetParticleFrameParameters(int i, double t,double t_step, Bunch* bunch) { - - + + double ta=2.418884326505e-17; //atomic unit of time double Ea=5.14220642e011; //Atomic unit of electric field double m=bunch->getMass(); @@ -257,14 +257,14 @@ part_t_step=t_step*coeff; //time step in frame of particle (in atomic units) -Ex_stat/=Ea; -Ey_stat/=Ea; +Ex_stat/=Ea; +Ey_stat/=Ea; Ez_stat/=Ea; - + } @@ -275,16 +275,16 @@ Ez_stat/=Ea; double HminusStripping::Gamma(double E_au) { - - + + double ta=2.418884326505e-17; //atomic unit of time double Ea=5.14220642e011; //Atomic unit of electric field - + double E = E_au*Ea*1e-8; double A = 7.96e-14; double C = 44.94; - - + + return ta/((A/E)*exp(C/E)); } @@ -293,32 +293,32 @@ double HminusStripping::Gamma(double E_au) { void HminusStripping::AmplSolver4step(double t_step, int i, Bunch* bunch) { - - + + double dt; - + if(bunch->flag(i) == 1) { - for(int j=1; j<5; j++) { - - if (j==4) dt=part_t_step; else dt=part_t_step/2.; - k_RungeKutt[j]=0; - k_RungeKutt[j] += (1 - (pop(i,0)+k_RungeKutt[j-1]*dt))*Gamma(Ez_stat); - + for(int j=1; j<5; j++) { + + if (j==4) dt=part_t_step; else dt=part_t_step/2.; + k_RungeKutt[j]=0; + k_RungeKutt[j] += (1 - (pop(i,0)+k_RungeKutt[j-1]*dt))*Gamma(Ez_stat); + } - - pop(i,0) += part_t_step*(k_RungeKutt[1]+2.*k_RungeKutt[2]+2.*k_RungeKutt[3]+k_RungeKutt[4])/6.; + + pop(i,0) += part_t_step*(k_RungeKutt[1]+2.*k_RungeKutt[2]+2.*k_RungeKutt[3]+k_RungeKutt[4])/6.; } - - + + if((index == 1)&&(pop(i,0) >= prob[i])) bunch->deleteParticleFast(i); - + } diff --git a/HminusStripping/HminusStripping.hh b/HminusStripping/HminusStripping.hh index 09bdcbc..df53648 100644 --- a/HminusStripping/HminusStripping.hh +++ b/HminusStripping/HminusStripping.hh @@ -5,9 +5,9 @@ #include "Python.h" -#include "ExternalEffects.hh" +#include #include "BaseLaserFieldSource.hh" -#include "ParticleAttributes.hh" +#include diff --git a/HminusStripping/wrap_hminus_stripping.cc b/HminusStripping/wrap_hminus_stripping.cc index e1a3e32..cbf2674 100644 --- a/HminusStripping/wrap_hminus_stripping.cc +++ b/HminusStripping/wrap_hminus_stripping.cc @@ -1,7 +1,6 @@ -#include "orbit_mpi.hh" -#include "pyORBIT_Object.hh" +#include +#include -#include "wrap_utils.hh" #include "wrap_hminus_stripping.hh" #include @@ -39,31 +38,31 @@ extern "C" { //initializator for python PyExternalEffects class //this is implementation of the __init__ method static int HminusStripping_init(pyORBIT_Object *self, PyObject *args, PyObject *kwds){ - - + + PyObject* pyStarkEffect=NULL; - + int method; int nVars = PyTuple_Size(args); - + if(!PyArg_ParseTuple(args,"i:",&method)) - {error("HminusStripping(method) - method = 0 or 1 param are needed");} - + {error("HminusStripping(method) - method = 0 or 1 param are needed");} + else { - + self->cpp_obj = new HminusStripping(method); ((HminusStripping*) self->cpp_obj)->setPyWrapper((PyObject*) self); } - - - + + + return 0; } - - - - // name([name]) - sets or returns the name of the External Effeects class + + + + // name([name]) - sets or returns the name of the External Effeects class static PyObject* HminusStripping_name(PyObject *self, PyObject *args){ HminusStripping* cpp_HminusStripping = (HminusStripping*) ((pyORBIT_Object*) self)->cpp_obj; const char* name = NULL; @@ -75,32 +74,32 @@ extern "C" { cpp_HminusStripping->setName(name_str); } return Py_BuildValue("s",cpp_HminusStripping->getName().c_str()); - } - - - - - - - - - - - - - - - - - - + } + + + + + + + + + + + + + + + + + + //----------------------------------------------------- //destructor for python PyExternalEffects class (__del__ method). //----------------------------------------------------- static void HminusStripping_del(pyORBIT_Object* self){ //std::cerr<<"The HminusStripping __del__ has been called!"<cpp_obj); - self->ob_type->tp_free((PyObject*)self); + Py_TYPE(self)->tp_free((PyObject*)self); } // defenition of the methods of the python PyExternalEffects wrapper class @@ -119,8 +118,7 @@ extern "C" { //new python PyExternalEffects wrapper type definition static PyTypeObject pyORBIT_HminusStripping_Type = { - PyObject_HEAD_INIT(NULL) - 0, /*ob_size*/ + PyVarObject_HEAD_INIT(NULL, 0) "HminusStripping", /*tp_name*/ sizeof(pyORBIT_Object), /*tp_basicsize*/ 0, /*tp_itemsize*/ @@ -158,10 +156,10 @@ extern "C" { (initproc) HminusStripping_init, /* tp_init */ 0, /* tp_alloc */ HminusStripping_new, /* tp_new */ - }; + }; + - //-------------------------------------------------- //Initialization function of the pyPyExternalEffects class //It will be called from Bunch wrapper initialization @@ -170,7 +168,7 @@ extern "C" { if (PyType_Ready(&pyORBIT_HminusStripping_Type) < 0) return; Py_INCREF(&pyORBIT_HminusStripping_Type); PyModule_AddObject(module, "HminusStripping", (PyObject *)&pyORBIT_HminusStripping_Type); - + } #ifdef __cplusplus diff --git a/HydrogenStarkParam/HydrogenStarkParam.cc b/HydrogenStarkParam/HydrogenStarkParam.cc index 57e7920..01b4bf1 100644 --- a/HydrogenStarkParam/HydrogenStarkParam.cc +++ b/HydrogenStarkParam/HydrogenStarkParam.cc @@ -7,19 +7,19 @@ // 04/21/2003 // // DESCRIPTION -// The base class for Python implementation of a field source. -// It should be sub-classed on Python level and implements +// The base class for Python implementation of a field source. +// It should be sub-classed on Python level and implements // getElectricField(x,y,z,t) and getMagneticField (x,y,z,t) methods. // The results of these methods will be available from the c++ level. // This is an example of embedding Python in C++ Orbit level. // // /////////////////////////////////////////////////////////////////////////// -#include "orbit_mpi.hh" -#include "BufferStore.hh" +#include +#include #include "HydrogenStarkParam.hh" -#include "LorentzTransformationEM.hh" +//#include "LorentzTransformationEM.hh" #include #include #include @@ -40,135 +40,135 @@ HydrogenStarkParam::HydrogenStarkParam(char* addressEG,int states) int rank_MPI,size_MPI; ORBIT_MPI_Comm_size(MPI_COMM_WORLD, &size_MPI); ORBIT_MPI_Comm_rank(MPI_COMM_WORLD, &rank_MPI); - + std::ifstream file; double F,alpha=7.297352570e-3; int k,ks,fi; char nameEG[MAX_LENGTH_ADDRESS]; - + st=states; - + levels=st*(1+st)*(1+2*st)/6; - - + + //this function measures parameter of input files (length, delta_F) using groung level file 000.txt - snprintf(nameEG,MAX_LENGTH_ADDRESS,"%s000.txt",addressEG); + snprintf(nameEG,MAX_LENGTH_ADDRESS,"%s000.txt",addressEG); if(rank_MPI == 0) { file.open(nameEG); file>>F>>F>>F>>delta_F; file.close(); file.open(nameEG);fi=0; while(!file.eof()) {file>>F>>F>>F;fi++;} file.close();n_data=fi-1; } - - ORBIT_MPI_Bcast(&delta_F, 1, MPI_DOUBLE,0,MPI_COMM_WORLD); + + ORBIT_MPI_Bcast(&delta_F, 1, MPI_DOUBLE,0,MPI_COMM_WORLD); ORBIT_MPI_Bcast(&n_data, 1, MPI_INT,0,MPI_COMM_WORLD); - + //allocating memory for dynamic massive of data that will be read fron files field_thresh = new double[levels+1]; - - energy=new double*[levels+1]; + + energy=new double*[levels+1]; for (int i=0;igetFreeDoubleArr(buff_index,n_data+10); - - + + //this loop reads energies and autoionization coefficients if(rank_MPI == 0) { for(int n=1;n>F>>energy[k][fi]>>gamma_autoionization[k][fi];} file_in.close(); - + file.open(nameEG); while(!file.eof()) file>>field_thresh[k]>>F>>F; file.close(); - + } } } } - - + + ORBIT_MPI_Bcast(field_thresh,levels+1 ,MPI_DOUBLE,0,MPI_COMM_WORLD); - - + + for(int n=1;n>F>>dipole_transition_x[k][ks][fi]>>dipole_transition_y[k][ks][fi]>>dipole_transition_z[k][ks][fi]; - + // this condition assumes that probability of spontaneous (see next line) and indused tansition between levels with the same principal quantum number n is sero if((ns==n)||(field_thresh[k]<=F)||(field_thresh[ks]<=F)) {dipole_transition_x[k][ks][fi]=0;dipole_transition_y[k][ks][fi]=0;dipole_transition_z[k][ks][fi]=0;} - + //this loop fills spontaneous relaxation (transition) of atom in relative atomic units gamma_spontaneous_relax[k][ks][fi]=fabs((4*alpha*alpha*alpha/3)*pow(energy[k][fi]-energy[ks][fi],3)*(pow(dipole_transition_x[k][ks][fi],2)+pow(dipole_transition_y[k][ks][fi],2)+pow(dipole_transition_z[k][ks][fi],2))); - + } file_in.close(); - + } } } @@ -180,59 +180,59 @@ HydrogenStarkParam::HydrogenStarkParam(char* addressEG,int states) for(int n=1;nsetUnusedDoubleArr(buff_index); + } + + BufferStore::getBufferStore()->setUnusedDoubleArr(buff_index); } HydrogenStarkParam::~HydrogenStarkParam() { - + for (int i=0;i1e-20)&&(gamma_autoionization[k][i+1]>1e-20)) Gamma=gamma_autoionization[k][i]*pow(gamma_autoionization[k][i+1]/gamma_autoionization[k][i],c); else Gamma=0; - + E=energy[k][i]+c*(energy[k][i+1]-energy[k][i]); - - + + } double HydrogenStarkParam::getStarkEnergy( - double mass, int n1, int n2, int m, + double mass, int n1, int n2, int m, double Ex,double Ey,double Ez, double Bx,double By,double Bz, double px,double py,double pz){ @@ -301,16 +301,16 @@ double Energy,Gamma; -LorentzTransformationEM::transform(mass,px,py,pz,Ex,Ey,Ez,Bx,By,Bz); +//LorentzTransformationEM::transform(mass,px,py,pz,Ex,Ey,Ez,Bx,By,Bz); -SetE(sqrt(Ex*Ex+Ey*Ey+Ez*Ez)/5.14220642e011); +SetE(sqrt(Ex*Ex+Ey*Ey+Ez*Ez)/5.14220642e011); GetEnergyAutoionization(convert3to1level(n1+n2+abs(m)+1,n1,m),Energy, Gamma); return Energy; } - - - - + + + + diff --git a/HydrogenStarkParam/HydrogenStarkParam.hh b/HydrogenStarkParam/HydrogenStarkParam.hh index a39e467..51d8b85 100644 --- a/HydrogenStarkParam/HydrogenStarkParam.hh +++ b/HydrogenStarkParam/HydrogenStarkParam.hh @@ -2,8 +2,8 @@ #define HYDROGENSTARKPARAM_HH_ -#include "tcomplex.hh" -#include "CppPyWrapper.hh" +#include +#include diff --git a/HydrogenStarkParam/wrap_hydrogen_stark_param.cc b/HydrogenStarkParam/wrap_hydrogen_stark_param.cc index efdb2b1..a0b46ec 100644 --- a/HydrogenStarkParam/wrap_hydrogen_stark_param.cc +++ b/HydrogenStarkParam/wrap_hydrogen_stark_param.cc @@ -1,7 +1,6 @@ -#include "orbit_mpi.hh" -#include "pyORBIT_Object.hh" +#include +#include -#include "wrap_utils.hh" #include "wrap_hydrogen_stark_param.hh" #include @@ -12,7 +11,6 @@ #include "HydrogenStarkParam.hh" //using namespace OrbitUtils; -using namespace wrap_orbit_utils; namespace wrap_hydrogen_stark_param{ @@ -45,26 +43,26 @@ extern "C" { if(!PyArg_ParseTuple( args,"si:",&addressEG,&states)){ error("HydrogenStarkParam(address,states) - params. are needed"); - } + } else { self->cpp_obj = new HydrogenStarkParam(addressEG,states); ((HydrogenStarkParam*) self->cpp_obj)->setPyWrapper((PyObject*) self); } - + return 0; - - + + } - - - - + + + + static PyObject* HydrogenStarkParam_getStarkEnergy(PyObject *self, PyObject *args){ HydrogenStarkParam* cpp_HydrogenStarkParam = (HydrogenStarkParam*)((pyORBIT_Object*) self)->cpp_obj; - - + + int nVars = PyTuple_Size(args); double val; double mass; @@ -80,7 +78,7 @@ extern "C" { double px; double py; double pz; - + //NO NEW OBJECT CREATED BY PyArg_ParseTuple! NO NEED OF Py_DECREF() if(!PyArg_ParseTuple( args,"diiiddddddddd:",&mass,&n1,&n2,&m,&E_x,&E_y,&E_z,&B_x,&B_y,&B_z,&px,&py,&pz)) @@ -89,12 +87,12 @@ extern "C" { val=cpp_HydrogenStarkParam->getStarkEnergy(mass,n1,n2,m,E_x,E_y,E_z,B_x,B_y,B_z,px,py,pz); return Py_BuildValue("d",val); } - - - - - + + + + + //----------------------------------------------------- //destructor for python PyBaseFieldSource class (__del__ method). @@ -102,7 +100,7 @@ extern "C" { static void HydrogenStarkParam_del(pyORBIT_Object* self){ //std::cerr<<"The CppBaseFieldSource __del__ has been called!"<cpp_obj); - self->ob_type->tp_free((PyObject*)self); + Py_TYPE(self)->tp_free((PyObject*)self); } // defenition of the methods of the python PyBaseFieldSource wrapper class @@ -128,8 +126,7 @@ extern "C" { //new python PyBaseFieldSource wrapper type definition static PyTypeObject pyORBIT_HydrogenStarkParam_Type = { - PyObject_HEAD_INIT(NULL) - 0, /*ob_size*/ + PyVarObject_HEAD_INIT(NULL, 0) "HydrogenStarkParam", /*tp_name*/ sizeof(pyORBIT_Object), /*tp_basicsize*/ 0, /*tp_itemsize*/ @@ -167,7 +164,7 @@ extern "C" { (initproc) HydrogenStarkParam_init, /* tp_init */ 0, /* tp_alloc */ HydrogenStarkParam_new, /* tp_new */ - }; + }; //-------------------------------------------------- //Initialization function of the pyPyBaseFieldSource class diff --git a/LS_utils/BaseLaserFieldSource.hh b/LS_utils/BaseLaserFieldSource.hh index 574e9fc..98e68c4 100644 --- a/LS_utils/BaseLaserFieldSource.hh +++ b/LS_utils/BaseLaserFieldSource.hh @@ -15,8 +15,8 @@ #ifndef BASELASERFIELDSOURCE_HH_ #define BASELASERFIELDSOURCE_HH_ -#include "CppPyWrapper.hh" -#include "tcomplex.hh" +#include +#include #define J tcomplex(0.,1.) diff --git a/LS_utils/FieldRotation.hh b/LS_utils/FieldRotation.hh index bdc23a0..66449b6 100644 --- a/LS_utils/FieldRotation.hh +++ b/LS_utils/FieldRotation.hh @@ -2,7 +2,7 @@ #define FIELDROTATION_HH_ -#include "tcomplex.hh" +#include class FieldRotation diff --git a/LS_utils/wrap_laserstripping.cc b/LS_utils/wrap_laserstripping.cc index 2f86494..0b127c2 100644 --- a/LS_utils/wrap_laserstripping.cc +++ b/LS_utils/wrap_laserstripping.cc @@ -1,4 +1,4 @@ -#include "orbit_mpi.hh" +#include #include "wrap_density_matrix.hh" #include "wrap_DM_noLaserField.hh" @@ -26,39 +26,39 @@ static PyMethodDef laserStrippingMethods[] = { {NULL,NULL} }; -#ifdef __cplusplus -extern "C" { -#endif +static struct PyModuleDef laserstrippingmodule = { + PyModuleDef_HEAD_INIT, + "laserstripping", + NULL, + -1, + laserStrippingMethods +}; - void initlaserstripping(){ - //create new module - PyObject* module = Py_InitModule("laserstripping",laserStrippingMethods); - wrap_density_matrix::initDensityMatrix(module); - wrap_DM_noLaserField::initDM_noLaserField(module); - wrap_schrodinger_equation::initSchrodingerEquation(module); - wrap_regular_grid_fs::initRegularGridFS(module); - wrap_hermite_gaussian_lf_mode::initHermiteGaussianLFmode(module); - wrap_las_field_container::initLaserFieldContainer(module); - wrap_hydrogen_stark_param::initHydrogenStarkParam(module); - wrap_two_level_atom::initTwoLevelAtom(module); - wrap_two_level_strong_field::initTwoLevelStrongField(module); - wrap_froissart_stora_lf::initFroissartStoraLF(module); - wrap_quad_em_field::initQuadEMfield(module); - wrap_const_em_field::initConstEMfield(module); - wrap_print_ext_effects::initPrintExtEffects(module); - wrap_record_evolution::initRecordEvolution(module); - wrap_walls::initWalls(module); - wrap_fringe_field::initFringeField(module); - wrap_stark::initStark(module); - wrap_stark_strong_field::initStarkStrongField(module); - wrap_cone_lf_mode::initConeLFmode(module); - wrap_hminus_stripping::initHminusStripping(module); - wrap_continuum_ss::initContinuumSS(module); - wrap_tdm_continuum::initTDMcontinuum(module); - wrap_snsstrippingmagnet::initSNSstrippingMagnet(module); - - } - -#ifdef __cplusplus +PyMODINIT_FUNC PyInit_laserstripping(void) { + PyObject* module = PyModule_Create(&laserstrippingmodule); + if (!module) return NULL; + wrap_density_matrix::initDensityMatrix(module); + wrap_DM_noLaserField::initDM_noLaserField(module); + wrap_schrodinger_equation::initSchrodingerEquation(module); + wrap_regular_grid_fs::initRegularGridFS(module); + wrap_hermite_gaussian_lf_mode::initHermiteGaussianLFmode(module); + wrap_las_field_container::initLaserFieldContainer(module); + wrap_hydrogen_stark_param::initHydrogenStarkParam(module); + wrap_two_level_atom::initTwoLevelAtom(module); + wrap_two_level_strong_field::initTwoLevelStrongField(module); + wrap_froissart_stora_lf::initFroissartStoraLF(module); + wrap_quad_em_field::initQuadEMfield(module); + wrap_const_em_field::initConstEMfield(module); + wrap_print_ext_effects::initPrintExtEffects(module); + wrap_record_evolution::initRecordEvolution(module); + wrap_walls::initWalls(module); + wrap_fringe_field::initFringeField(module); + wrap_stark::initStark(module); + wrap_stark_strong_field::initStarkStrongField(module); + wrap_cone_lf_mode::initConeLFmode(module); + wrap_hminus_stripping::initHminusStripping(module); + wrap_continuum_ss::initContinuumSS(module); + wrap_tdm_continuum::initTDMcontinuum(module); + wrap_snsstrippingmagnet::initSNSstrippingMagnet(module); + return module; } -#endif diff --git a/LaserFieldContainer/LaserFieldContainer.cc b/LaserFieldContainer/LaserFieldContainer.cc index 9868bf2..a7ba348 100644 --- a/LaserFieldContainer/LaserFieldContainer.cc +++ b/LaserFieldContainer/LaserFieldContainer.cc @@ -7,8 +7,8 @@ // 04/21/2003 // // DESCRIPTION -// The base class for Python implementation of a field source. -// It should be sub-classed on Python level and implements +// The base class for Python implementation of a field source. +// It should be sub-classed on Python level and implements // getElectricField(x,y,z,t) and getMagneticField (x,y,z,t) methods. // The results of these methods will be available from the c++ level. // This is an example of embedding Python in C++ Orbit level. @@ -17,8 +17,8 @@ /////////////////////////////////////////////////////////////////////////// #include "LaserFieldContainer.hh" -#include "orbit_mpi.hh" -#include "OrbitConst.hh" +#include +#include #include @@ -40,16 +40,16 @@ LaserFieldContainer::~LaserFieldContainer() void LaserFieldContainer::getLaserElectricField(double x, double y, double z, double t, tcomplex& E_x, tcomplex& E_y, tcomplex& E_z){ - + } void LaserFieldContainer::getLaserMagneticField(double x, double y, double z, double t, tcomplex& B_x, tcomplex& B_y, tcomplex& B_z){ - + } double LaserFieldContainer::getFrequencyOmega(double m,double x, double y, double z, double px, double py, double pz, double t){ - + return 0; } diff --git a/LaserFieldContainer/wrap_las_field_container.cc b/LaserFieldContainer/wrap_las_field_container.cc index 9003186..8906906 100644 --- a/LaserFieldContainer/wrap_las_field_container.cc +++ b/LaserFieldContainer/wrap_las_field_container.cc @@ -1,7 +1,6 @@ -#include "orbit_mpi.hh" -#include "pyORBIT_Object.hh" +#include +#include -#include "wrap_utils.hh" #include "wrap_las_field_container.hh" #include @@ -10,7 +9,6 @@ #include "LaserFieldContainer.hh" using namespace OrbitUtils; -using namespace wrap_orbit_utils; namespace wrap_las_field_container{ @@ -40,11 +38,11 @@ extern "C" { self->cpp_obj = new LaserFieldContainer(); //constructor with parameters is nedeed return 0; } - - - - + + + + //----------------------------------------------------- //destructor for python PyBaseFieldSource class (__del__ method). @@ -52,13 +50,13 @@ extern "C" { static void LaserFieldContainer_del(pyORBIT_Object* self){ //std::cerr<<"The CppBaseFieldSource __del__ has been called!"<cpp_obj); - self->ob_type->tp_free((PyObject*)self); + Py_TYPE(self)->tp_free((PyObject*)self); } // defenition of the methods of the python PyBaseFieldSource wrapper class // they will be vailable from python level static PyMethodDef LaserFieldContainerClassMethods[] = { - + {NULL} }; @@ -70,8 +68,7 @@ extern "C" { //new python PyBaseFieldSource wrapper type definition static PyTypeObject pyORBIT_LaserFieldContainer_Type = { - PyObject_HEAD_INIT(NULL) - 0, /*ob_size*/ + PyVarObject_HEAD_INIT(NULL, 0) "laserFieldContainer", /*tp_name*/ sizeof(pyORBIT_Object), /*tp_basicsize*/ 0, /*tp_itemsize*/ @@ -109,7 +106,7 @@ extern "C" { (initproc) LaserFieldContainer_init, /* tp_init */ 0, /* tp_alloc */ LaserFieldContainer_new, /* tp_new */ - }; + }; //-------------------------------------------------- //Initialization function of the pyPyBaseFieldSource class diff --git a/Quad_EM_field/ConstEMfield.cc b/Quad_EM_field/ConstEMfield.cc index 4f52de7..ffa8a35 100644 --- a/Quad_EM_field/ConstEMfield.cc +++ b/Quad_EM_field/ConstEMfield.cc @@ -7,8 +7,8 @@ // 04/21/2003 // // DESCRIPTION -// The base class for Python implementation of a field source. -// It should be sub-classed on Python level and implements +// The base class for Python implementation of a field source. +// It should be sub-classed on Python level and implements // getElectricField(x,y,z,t) and getMagneticField (x,y,z,t) methods. // The results of these methods will be available from the c++ level. // This is an example of embedding Python in C++ Orbit level. @@ -17,35 +17,35 @@ /////////////////////////////////////////////////////////////////////////// #include "ConstEMfield.hh" -#include "orbit_mpi.hh" +#include #include using namespace OrbitUtils; ConstEMfield::ConstEMfield() { - + ExConst=0; EyConst=0; EzConst=0; - + BxConst=0; ByConst=0; BzConst=0; - + } ConstEMfield::ConstEMfield(double E_x,double E_y,double E_z,double B_x,double B_y,double B_z) { - + ExConst=E_x; EyConst=E_y; EzConst=E_z; - + BxConst=B_x; ByConst=B_y; BzConst=B_z; - + } @@ -55,14 +55,14 @@ ConstEMfield::~ConstEMfield() } -void ConstEMfield::getElectricMagneticField(double x, double y, double z, double t, +void ConstEMfield::getElectricMagneticField(double x, double y, double z, double t, double& E_x, double& E_y, double& E_z, double& B_x, double& B_y, double& B_z) { - - E_x = ExConst; E_y = EyConst; E_z = EzConst; - B_x = BxConst; B_y = ByConst; B_z = BzConst; + + E_x = ExConst; E_y = EyConst; E_z = EzConst; + B_x = BxConst; B_y = ByConst; B_z = BzConst; } diff --git a/Quad_EM_field/ConstEMfield.hh b/Quad_EM_field/ConstEMfield.hh index 107b944..1e27fc8 100644 --- a/Quad_EM_field/ConstEMfield.hh +++ b/Quad_EM_field/ConstEMfield.hh @@ -24,7 +24,7 @@ #include "Python.h" -#include "BaseFieldSource.hh" +#include namespace OrbitUtils{ diff --git a/Quad_EM_field/QuadEMfield.cc b/Quad_EM_field/QuadEMfield.cc index ece67ee..6bb26c5 100644 --- a/Quad_EM_field/QuadEMfield.cc +++ b/Quad_EM_field/QuadEMfield.cc @@ -7,8 +7,8 @@ // 04/21/2003 // // DESCRIPTION -// The base class for Python implementation of a field source. -// It should be sub-classed on Python level and implements +// The base class for Python implementation of a field source. +// It should be sub-classed on Python level and implements // getElectricField(x,y,z,t) and getMagneticField (x,y,z,t) methods. // The results of these methods will be available from the c++ level. // This is an example of embedding Python in C++ Orbit level. @@ -17,7 +17,7 @@ /////////////////////////////////////////////////////////////////////////// #include "QuadEMfield.hh" -#include "orbit_mpi.hh" +#include #include using namespace OrbitUtils; @@ -27,35 +27,35 @@ QuadEMfield::QuadEMfield() Ex0v=0; Ey0v=0; Ez0v=0; - + Bx0v=0; By0v=0; Bz0v=0; - + x0v=0; y0v=0; z0v=0; - + cxExv=0; cyExv=0; czExv=0; - + cxEyv=0; cyEyv=0; czEyv=0; - + cxEzv=0; cyEzv=0; czEzv=0; - + cxBxv=0; cyBxv=0; czBxv=0; - + cxByv=0; cyByv=0; czByv=0; - + cxBzv=0; cyBzv=0; czBzv=0; @@ -68,26 +68,26 @@ QuadEMfield::~QuadEMfield() -void QuadEMfield::getElectricMagneticField(double x, double y, double z, double t, +void QuadEMfield::getElectricMagneticField(double x, double y, double z, double t, double& E_x, double& E_y, double& E_z, double& B_x, double& B_y, double& B_z) -{ +{ x-=x0v; y-=y0v; z-=z0v; - - E_x = Ex0v+cxExv*x+cyExv*y+czExv*z; - E_y = Ey0v+cxEyv*x+cyEyv*y+czEyv*z; - E_z = Ez0v+cxEzv*x+cyEzv*y+czEzv*z; - - B_x = Bx0v+cxBxv*x+cyBxv*y+czBxv*z; - B_y = By0v+cxByv*x+cyByv*y+czByv*z; + + E_x = Ex0v+cxExv*x+cyExv*y+czExv*z; + E_y = Ey0v+cxEyv*x+cyEyv*y+czEyv*z; + E_z = Ez0v+cxEzv*x+cyEzv*y+czEzv*z; + + B_x = Bx0v+cxBxv*x+cyBxv*y+czBxv*z; + B_y = By0v+cxByv*x+cyByv*y+czByv*z; B_z = Bz0v+cxBzv*x+cyBzv*y+czBzv*z; - + return; } diff --git a/Quad_EM_field/QuadEMfield.hh b/Quad_EM_field/QuadEMfield.hh index e81c8b5..a68b9be 100644 --- a/Quad_EM_field/QuadEMfield.hh +++ b/Quad_EM_field/QuadEMfield.hh @@ -24,7 +24,7 @@ #include "Python.h" -#include "BaseFieldSource.hh" +#include namespace OrbitUtils{ diff --git a/Quad_EM_field/wrap_const_em_field.cc b/Quad_EM_field/wrap_const_em_field.cc index c8c35bc..9161b15 100644 --- a/Quad_EM_field/wrap_const_em_field.cc +++ b/Quad_EM_field/wrap_const_em_field.cc @@ -1,7 +1,6 @@ -#include "orbit_mpi.hh" -#include "pyORBIT_Object.hh" +#include +#include -#include "wrap_utils.hh" #include "wrap_const_em_field.hh" #include @@ -10,7 +9,6 @@ #include "ConstEMfield.hh" using namespace OrbitUtils; -using namespace wrap_orbit_utils; namespace wrap_const_em_field{ @@ -37,111 +35,111 @@ extern "C" { //initializator for python CppBaseFieldSource class //this is implementation of the __init__ method static int ConstEMfield_init(pyORBIT_Object *self, PyObject *args, PyObject *kwds){ - + // self->cpp_obj = new LSFieldSource(); - + double E_x; double E_y; double E_z; double B_x; double B_y; double B_z; - + int nVars = PyTuple_Size(args); - + // if (PyArg_ParseTuple( args,"")) {self->cpp_obj = new LSFieldSource();} // if(PyArg_ParseTuple(args,"dddddd:",&E_x,&E_y,&E_z,&B_x,&B_y,&B_z)) { self->cpp_obj = new LSFieldSource(E_x,E_y,E_z,B_x,B_y,B_z);} - + if(nVars==0) if (!PyArg_ParseTuple( args,"")) {} else {self->cpp_obj = new ConstEMfield();} - if(nVars==6) if (!PyArg_ParseTuple(args,"dddddd:",&E_x,&E_y,&E_z,&B_x,&B_y,&B_z)) + if(nVars==6) if (!PyArg_ParseTuple(args,"dddddd:",&E_x,&E_y,&E_z,&B_x,&B_y,&B_z)) {error("Parameters (E_x,E_y,E_z,B_x,B_y,B_z) -are needed");} else {self->cpp_obj = new ConstEMfield(E_x,E_y,E_z,B_x,B_y,B_z);} ((ConstEMfield*) self->cpp_obj)->setPyWrapper((PyObject*) self); - + return 0; } - - + + static PyObject* ConstEMfield_Ex(PyObject *self, PyObject *args){ - ConstEMfield* cpp_ConstEMfield = (ConstEMfield*)((pyORBIT_Object*) self)->cpp_obj; + ConstEMfield* cpp_ConstEMfield = (ConstEMfield*)((pyORBIT_Object*) self)->cpp_obj; double Ex; //NO NEW OBJECT CREATED BY PyArg_ParseTuple! NO NEED OF Py_DECREF() if(!PyArg_ParseTuple( args,"d:",&Ex)) error(" SetupPrint(Ex - parameter is needed"); - else + else cpp_ConstEMfield->Ex(Ex); - + Py_INCREF(Py_None); return Py_None; } - + static PyObject* ConstEMfield_Ey(PyObject *self, PyObject *args){ - ConstEMfield* cpp_ConstEMfield = (ConstEMfield*)((pyORBIT_Object*) self)->cpp_obj; + ConstEMfield* cpp_ConstEMfield = (ConstEMfield*)((pyORBIT_Object*) self)->cpp_obj; double Ey; //NO NEW OBJECT CREATED BY PyArg_ParseTuple! NO NEED OF Py_DECREF() if(!PyArg_ParseTuple( args,"d:",&Ey)) error(" SetupPrint(Ey - parameter is needed"); - else + else cpp_ConstEMfield->Ey(Ey); - + Py_INCREF(Py_None); return Py_None; } - + static PyObject* ConstEMfield_Ez(PyObject *self, PyObject *args){ - ConstEMfield* cpp_ConstEMfield = (ConstEMfield*)((pyORBIT_Object*) self)->cpp_obj; + ConstEMfield* cpp_ConstEMfield = (ConstEMfield*)((pyORBIT_Object*) self)->cpp_obj; double Ez; //NO NEW OBJECT CREATED BY PyArg_ParseTuple! NO NEED OF Py_DECREF() if(!PyArg_ParseTuple( args,"d:",&Ez)) error(" SetupPrint(Ez - parameter is needed"); - else + else cpp_ConstEMfield->Ez(Ez); - + Py_INCREF(Py_None); return Py_None; } - + static PyObject* ConstEMfield_Bx(PyObject *self, PyObject *args){ - ConstEMfield* cpp_ConstEMfield = (ConstEMfield*)((pyORBIT_Object*) self)->cpp_obj; + ConstEMfield* cpp_ConstEMfield = (ConstEMfield*)((pyORBIT_Object*) self)->cpp_obj; double Bx; //NO NEW OBJECT CREATED BY PyArg_ParseTuple! NO NEED OF Py_DECREF() if(!PyArg_ParseTuple( args,"d:",&Bx)) error(" SetupPrint(Bx - parameter is needed"); - else + else cpp_ConstEMfield->Bx(Bx); - + Py_INCREF(Py_None); return Py_None; } - + static PyObject* ConstEMfield_By(PyObject *self, PyObject *args){ - ConstEMfield* cpp_ConstEMfield = (ConstEMfield*)((pyORBIT_Object*) self)->cpp_obj; + ConstEMfield* cpp_ConstEMfield = (ConstEMfield*)((pyORBIT_Object*) self)->cpp_obj; double By; //NO NEW OBJECT CREATED BY PyArg_ParseTuple! NO NEED OF Py_DECREF() if(!PyArg_ParseTuple( args,"d:",&By)) error(" SetupPrint(By - parameter is needed"); - else + else cpp_ConstEMfield->By(By); - + Py_INCREF(Py_None); return Py_None; } - + static PyObject* ConstEMfield_Bz(PyObject *self, PyObject *args){ - ConstEMfield* cpp_ConstEMfield = (ConstEMfield*)((pyORBIT_Object*) self)->cpp_obj; + ConstEMfield* cpp_ConstEMfield = (ConstEMfield*)((pyORBIT_Object*) self)->cpp_obj; double Bz; //NO NEW OBJECT CREATED BY PyArg_ParseTuple! NO NEED OF Py_DECREF() if(!PyArg_ParseTuple( args,"d:",&Bz)) error(" SetupPrint(Bz - parameter is needed"); - else + else cpp_ConstEMfield->Bz(Bz); - + Py_INCREF(Py_None); return Py_None; } - - + + //----------------------------------------------------- @@ -150,7 +148,7 @@ extern "C" { static void ConstEMfield_del(pyORBIT_Object* self){ //std::cerr<<"The ConstEMfield __del__ has been called!"<cpp_obj); - self->ob_type->tp_free((PyObject*)self); + Py_TYPE(self)->tp_free((PyObject*)self); } // defenition of the methods of the python PyBaseFieldSource wrapper class @@ -162,7 +160,7 @@ extern "C" { { "Bx", ConstEMfield_Bx, METH_VARARGS,"Sets parameter Bx of EM field."}, { "By", ConstEMfield_By, METH_VARARGS,"Sets parameter By of EM field."}, { "Bz", ConstEMfield_Bz, METH_VARARGS,"Sets parameter Bz of EM field."}, - + {NULL} }; @@ -174,8 +172,7 @@ extern "C" { //new python PyBaseFieldSource wrapper type definition static PyTypeObject pyORBIT_ConstEMfield_Type = { - PyObject_HEAD_INIT(NULL) - 0, /*ob_size*/ + PyVarObject_HEAD_INIT(NULL, 0) "ConstEMfield", /*tp_name*/ sizeof(pyORBIT_Object), /*tp_basicsize*/ 0, /*tp_itemsize*/ @@ -213,7 +210,7 @@ extern "C" { (initproc) ConstEMfield_init, /* tp_init */ 0, /* tp_alloc */ ConstEMfield_new, /* tp_new */ - }; + }; //-------------------------------------------------- //Initialization function of the pyPyBaseFieldSource class diff --git a/Quad_EM_field/wrap_quad_em_field.cc b/Quad_EM_field/wrap_quad_em_field.cc index dc18a89..1a97950 100644 --- a/Quad_EM_field/wrap_quad_em_field.cc +++ b/Quad_EM_field/wrap_quad_em_field.cc @@ -1,7 +1,6 @@ -#include "orbit_mpi.hh" -#include "pyORBIT_Object.hh" +#include +#include -#include "wrap_utils.hh" #include "wrap_quad_em_field.hh" #include @@ -10,7 +9,6 @@ #include "QuadEMfield.hh" using namespace OrbitUtils; -using namespace wrap_orbit_utils; namespace wrap_quad_em_field{ @@ -37,377 +35,377 @@ extern "C" { //initializator for python CppBaseFieldSource class //this is implementation of the __init__ method static int QuadEMfield_init(pyORBIT_Object *self, PyObject *args, PyObject *kwds){ - + self->cpp_obj = new QuadEMfield(); - + ((QuadEMfield*) self->cpp_obj)->setPyWrapper((PyObject*) self); - + return 0; } - - + + static PyObject* QuadEMfield_x0(PyObject *self, PyObject *args){ - QuadEMfield* cpp_QuadEMfield = (QuadEMfield*)((pyORBIT_Object*) self)->cpp_obj; + QuadEMfield* cpp_QuadEMfield = (QuadEMfield*)((pyORBIT_Object*) self)->cpp_obj; double x0; //NO NEW OBJECT CREATED By0 PyArg_ParseTuple! NO NEED OF Py_DECREF() if(!PyArg_ParseTuple( args,"d:",&x0)) error(" SetupPrint(x0 - parameter is needed"); - else + else cpp_QuadEMfield->x0(x0); - + Py_INCREF(Py_None); return Py_None; } - + static PyObject* QuadEMfield_y0(PyObject *self, PyObject *args){ - QuadEMfield* cpp_QuadEMfield = (QuadEMfield*)((pyORBIT_Object*) self)->cpp_obj; + QuadEMfield* cpp_QuadEMfield = (QuadEMfield*)((pyORBIT_Object*) self)->cpp_obj; double y0; //NO NEW OBJECT CREATED By0 PyArg_ParseTuple! NO NEED OF Py_DECREF() if(!PyArg_ParseTuple( args,"d:",&y0)) error(" SetupPrint(y0 - parameter is needed"); - else + else cpp_QuadEMfield->y0(y0); - + Py_INCREF(Py_None); return Py_None; } - + static PyObject* QuadEMfield_z0(PyObject *self, PyObject *args){ - QuadEMfield* cpp_QuadEMfield = (QuadEMfield*)((pyORBIT_Object*) self)->cpp_obj; + QuadEMfield* cpp_QuadEMfield = (QuadEMfield*)((pyORBIT_Object*) self)->cpp_obj; double z0; //NO NEW OBJECT CREATED By0 PyArg_ParseTuple! NO NEED OF Py_DECREF() if(!PyArg_ParseTuple( args,"d:",&z0)) error(" SetupPrint(z0 - parameter is needed"); - else + else cpp_QuadEMfield->z0(z0); - + Py_INCREF(Py_None); return Py_None; } - - - + + + static PyObject* QuadEMfield_Ex0(PyObject *self, PyObject *args){ - QuadEMfield* cpp_QuadEMfield = (QuadEMfield*)((pyORBIT_Object*) self)->cpp_obj; + QuadEMfield* cpp_QuadEMfield = (QuadEMfield*)((pyORBIT_Object*) self)->cpp_obj; double Ex0; //NO NEW OBJECT CREATED By0 PyArg_ParseTuple! NO NEED OF Py_DECREF() if(!PyArg_ParseTuple( args,"d:",&Ex0)) error(" SetupPrint(Ex0 - parameter is needed"); - else + else cpp_QuadEMfield->Ex0(Ex0); - + Py_INCREF(Py_None); return Py_None; } - + static PyObject* QuadEMfield_Ey0(PyObject *self, PyObject *args){ - QuadEMfield* cpp_QuadEMfield = (QuadEMfield*)((pyORBIT_Object*) self)->cpp_obj; + QuadEMfield* cpp_QuadEMfield = (QuadEMfield*)((pyORBIT_Object*) self)->cpp_obj; double Ey0; //NO NEW OBJECT CREATED By0 PyArg_ParseTuple! NO NEED OF Py_DECREF() if(!PyArg_ParseTuple( args,"d:",&Ey0)) error(" SetupPrint(Ey0 - parameter is needed"); - else + else cpp_QuadEMfield->Ey0(Ey0); - + Py_INCREF(Py_None); return Py_None; } - + static PyObject* QuadEMfield_Ez0(PyObject *self, PyObject *args){ - QuadEMfield* cpp_QuadEMfield = (QuadEMfield*)((pyORBIT_Object*) self)->cpp_obj; + QuadEMfield* cpp_QuadEMfield = (QuadEMfield*)((pyORBIT_Object*) self)->cpp_obj; double Ez0; //NO NEW OBJECT CREATED By0 PyArg_ParseTuple! NO NEED OF Py_DECREF() if(!PyArg_ParseTuple( args,"d:",&Ez0)) error(" SetupPrint(Ez0 - parameter is needed"); - else + else cpp_QuadEMfield->Ez0(Ez0); - + Py_INCREF(Py_None); return Py_None; } - + static PyObject* QuadEMfield_Bx0(PyObject *self, PyObject *args){ - QuadEMfield* cpp_QuadEMfield = (QuadEMfield*)((pyORBIT_Object*) self)->cpp_obj; + QuadEMfield* cpp_QuadEMfield = (QuadEMfield*)((pyORBIT_Object*) self)->cpp_obj; double Bx0; //NO NEW OBJECT CREATED By0 PyArg_ParseTuple! NO NEED OF Py_DECREF() if(!PyArg_ParseTuple( args,"d:",&Bx0)) error(" SetupPrint(Bx0 - parameter is needed"); - else + else cpp_QuadEMfield->Bx0(Bx0); - + Py_INCREF(Py_None); return Py_None; } - + static PyObject* QuadEMfield_By0(PyObject *self, PyObject *args){ - QuadEMfield* cpp_QuadEMfield = (QuadEMfield*)((pyORBIT_Object*) self)->cpp_obj; + QuadEMfield* cpp_QuadEMfield = (QuadEMfield*)((pyORBIT_Object*) self)->cpp_obj; double By0; //NO NEW OBJECT CREATED By0 PyArg_ParseTuple! NO NEED OF Py_DECREF() if(!PyArg_ParseTuple( args,"d:",&By0)) error(" SetupPrint(By0 - parameter is needed"); - else + else cpp_QuadEMfield->By0(By0); - + Py_INCREF(Py_None); return Py_None; } - + static PyObject* QuadEMfield_Bz0(PyObject *self, PyObject *args){ - QuadEMfield* cpp_QuadEMfield = (QuadEMfield*)((pyORBIT_Object*) self)->cpp_obj; + QuadEMfield* cpp_QuadEMfield = (QuadEMfield*)((pyORBIT_Object*) self)->cpp_obj; double Bz0; //NO NEW OBJECT CREATED By0 PyArg_ParseTuple! NO NEED OF Py_DECREF() if(!PyArg_ParseTuple( args,"d:",&Bz0)) error(" SetupPrint(Bz0 - parameter is needed"); - else + else cpp_QuadEMfield->Bz0(Bz0); - + Py_INCREF(Py_None); return Py_None; } - + static PyObject* QuadEMfield_cxEx(PyObject *self, PyObject *args){ - QuadEMfield* cpp_QuadEMfield = (QuadEMfield*)((pyORBIT_Object*) self)->cpp_obj; + QuadEMfield* cpp_QuadEMfield = (QuadEMfield*)((pyORBIT_Object*) self)->cpp_obj; double cxEx; //NO NEW OBJECT CREATED By0 PyArg_ParseTuple! NO NEED OF Py_DECREF() if(!PyArg_ParseTuple( args,"d:",&cxEx)) error(" SetupPrint(Bz0 - parameter is needed"); - else + else cpp_QuadEMfield->cxEx(cxEx); - + Py_INCREF(Py_None); return Py_None; } - + static PyObject* QuadEMfield_cyEx(PyObject *self, PyObject *args){ - QuadEMfield* cpp_QuadEMfield = (QuadEMfield*)((pyORBIT_Object*) self)->cpp_obj; + QuadEMfield* cpp_QuadEMfield = (QuadEMfield*)((pyORBIT_Object*) self)->cpp_obj; double cyEx; //NO NEW OBJECT CREATED By0 PyArg_ParseTuple! NO NEED OF Py_DECREF() if(!PyArg_ParseTuple( args,"d:",&cyEx)) error(" SetupPrint(Bz0 - parameter is needed"); - else + else cpp_QuadEMfield->cyEx(cyEx); - + Py_INCREF(Py_None); return Py_None; } - + static PyObject* QuadEMfield_czEx(PyObject *self, PyObject *args){ - QuadEMfield* cpp_QuadEMfield = (QuadEMfield*)((pyORBIT_Object*) self)->cpp_obj; + QuadEMfield* cpp_QuadEMfield = (QuadEMfield*)((pyORBIT_Object*) self)->cpp_obj; double czEx; //NO NEW OBJECT CREATED By0 PyArg_ParseTuple! NO NEED OF Py_DECREF() if(!PyArg_ParseTuple( args,"d:",&czEx)) error(" SetupPrint(Bz0 - parameter is needed"); - else + else cpp_QuadEMfield->czEx(czEx); - + Py_INCREF(Py_None); return Py_None; } - + static PyObject* QuadEMfield_cxEy(PyObject *self, PyObject *args){ - QuadEMfield* cpp_QuadEMfield = (QuadEMfield*)((pyORBIT_Object*) self)->cpp_obj; + QuadEMfield* cpp_QuadEMfield = (QuadEMfield*)((pyORBIT_Object*) self)->cpp_obj; double cxEy; //NO NEW OBJECT CREATED By0 PyArg_ParseTuple! NO NEED OF Py_DECREF() if(!PyArg_ParseTuple( args,"d:",&cxEy)) error(" SetupPrint(Bz0 - parameter is needed"); - else + else cpp_QuadEMfield->cxEy(cxEy); - + Py_INCREF(Py_None); return Py_None; } - + static PyObject* QuadEMfield_cyEy(PyObject *self, PyObject *args){ - QuadEMfield* cpp_QuadEMfield = (QuadEMfield*)((pyORBIT_Object*) self)->cpp_obj; + QuadEMfield* cpp_QuadEMfield = (QuadEMfield*)((pyORBIT_Object*) self)->cpp_obj; double cyEy; //NO NEW OBJECT CREATED By0 PyArg_ParseTuple! NO NEED OF Py_DECREF() if(!PyArg_ParseTuple( args,"d:",&cyEy)) error(" SetupPrint(Bz0 - parameter is needed"); - else + else cpp_QuadEMfield->cyEy(cyEy); - + Py_INCREF(Py_None); return Py_None; } - + static PyObject* QuadEMfield_czEy(PyObject *self, PyObject *args){ - QuadEMfield* cpp_QuadEMfield = (QuadEMfield*)((pyORBIT_Object*) self)->cpp_obj; + QuadEMfield* cpp_QuadEMfield = (QuadEMfield*)((pyORBIT_Object*) self)->cpp_obj; double czEy; //NO NEW OBJECT CREATED By0 PyArg_ParseTuple! NO NEED OF Py_DECREF() if(!PyArg_ParseTuple( args,"d:",&czEy)) error(" SetupPrint(Bz0 - parameter is needed"); - else + else cpp_QuadEMfield->czEy(czEy); - + Py_INCREF(Py_None); return Py_None; } - + static PyObject* QuadEMfield_cxEz(PyObject *self, PyObject *args){ - QuadEMfield* cpp_QuadEMfield = (QuadEMfield*)((pyORBIT_Object*) self)->cpp_obj; + QuadEMfield* cpp_QuadEMfield = (QuadEMfield*)((pyORBIT_Object*) self)->cpp_obj; double cxEz; //NO NEW OBJECT CREATED By0 PyArg_ParseTuple! NO NEED OF Py_DECREF() if(!PyArg_ParseTuple( args,"d:",&cxEz)) error(" SetupPrint(Bz0 - parameter is needed"); - else + else cpp_QuadEMfield->cxEz(cxEz); - + Py_INCREF(Py_None); return Py_None; } - + static PyObject* QuadEMfield_cyEz(PyObject *self, PyObject *args){ - QuadEMfield* cpp_QuadEMfield = (QuadEMfield*)((pyORBIT_Object*) self)->cpp_obj; + QuadEMfield* cpp_QuadEMfield = (QuadEMfield*)((pyORBIT_Object*) self)->cpp_obj; double cyEz; //NO NEW OBJECT CREATED By0 PyArg_ParseTuple! NO NEED OF Py_DECREF() if(!PyArg_ParseTuple( args,"d:",&cyEz)) error(" SetupPrint(Bz0 - parameter is needed"); - else + else cpp_QuadEMfield->cyEz(cyEz); - + Py_INCREF(Py_None); return Py_None; } - + static PyObject* QuadEMfield_czEz(PyObject *self, PyObject *args){ - QuadEMfield* cpp_QuadEMfield = (QuadEMfield*)((pyORBIT_Object*) self)->cpp_obj; + QuadEMfield* cpp_QuadEMfield = (QuadEMfield*)((pyORBIT_Object*) self)->cpp_obj; double czEz; //NO NEW OBJECT CREATED By0 PyArg_ParseTuple! NO NEED OF Py_DECREF() if(!PyArg_ParseTuple( args,"d:",&czEz)) error(" SetupPrint(Bz0 - parameter is needed"); - else + else cpp_QuadEMfield->czEz(czEz); - + Py_INCREF(Py_None); return Py_None; } - + static PyObject* QuadEMfield_cxBx(PyObject *self, PyObject *args){ - QuadEMfield* cpp_QuadEMfield = (QuadEMfield*)((pyORBIT_Object*) self)->cpp_obj; + QuadEMfield* cpp_QuadEMfield = (QuadEMfield*)((pyORBIT_Object*) self)->cpp_obj; double cxBx; //NO NEW OBJECT CREATED By0 PyArg_ParseTuple! NO NEED OF Py_DECREF() if(!PyArg_ParseTuple( args,"d:",&cxBx)) error(" SetupPrint(Bz0 - parameter is needed"); - else + else cpp_QuadEMfield->cxBx(cxBx); - + Py_INCREF(Py_None); return Py_None; } - + static PyObject* QuadEMfield_cyBx(PyObject *self, PyObject *args){ - QuadEMfield* cpp_QuadEMfield = (QuadEMfield*)((pyORBIT_Object*) self)->cpp_obj; + QuadEMfield* cpp_QuadEMfield = (QuadEMfield*)((pyORBIT_Object*) self)->cpp_obj; double cyBx; //NO NEW OBJECT CREATED By0 PyArg_ParseTuple! NO NEED OF Py_DECREF() if(!PyArg_ParseTuple( args,"d:",&cyBx)) error(" SetupPrint(Bz0 - parameter is needed"); - else + else cpp_QuadEMfield->cyBx(cyBx); - + Py_INCREF(Py_None); return Py_None; } - + static PyObject* QuadEMfield_czBx(PyObject *self, PyObject *args){ - QuadEMfield* cpp_QuadEMfield = (QuadEMfield*)((pyORBIT_Object*) self)->cpp_obj; + QuadEMfield* cpp_QuadEMfield = (QuadEMfield*)((pyORBIT_Object*) self)->cpp_obj; double czBx; //NO NEW OBJECT CREATED By0 PyArg_ParseTuple! NO NEED OF Py_DECREF() if(!PyArg_ParseTuple( args,"d:",&czBx)) error(" SetupPrint(Bz0 - parameter is needed"); - else + else cpp_QuadEMfield->czBx(czBx); - + Py_INCREF(Py_None); return Py_None; } - + static PyObject* QuadEMfield_cxBy(PyObject *self, PyObject *args){ - QuadEMfield* cpp_QuadEMfield = (QuadEMfield*)((pyORBIT_Object*) self)->cpp_obj; + QuadEMfield* cpp_QuadEMfield = (QuadEMfield*)((pyORBIT_Object*) self)->cpp_obj; double cxBy; //NO NEW OBJECT CREATED By0 PyArg_ParseTuple! NO NEED OF Py_DECREF() if(!PyArg_ParseTuple( args,"d:",&cxBy)) error(" SetupPrint(Bz0 - parameter is needed"); - else + else cpp_QuadEMfield->cxBy(cxBy); - + Py_INCREF(Py_None); return Py_None; } - + static PyObject* QuadEMfield_cyBy(PyObject *self, PyObject *args){ - QuadEMfield* cpp_QuadEMfield = (QuadEMfield*)((pyORBIT_Object*) self)->cpp_obj; + QuadEMfield* cpp_QuadEMfield = (QuadEMfield*)((pyORBIT_Object*) self)->cpp_obj; double cyBy; //NO NEW OBJECT CREATED By0 PyArg_ParseTuple! NO NEED OF Py_DECREF() if(!PyArg_ParseTuple( args,"d:",&cyBy)) error(" SetupPrint(Bz0 - parameter is needed"); - else + else cpp_QuadEMfield->cyBy(cyBy); - + Py_INCREF(Py_None); return Py_None; } - + static PyObject* QuadEMfield_czBy(PyObject *self, PyObject *args){ - QuadEMfield* cpp_QuadEMfield = (QuadEMfield*)((pyORBIT_Object*) self)->cpp_obj; + QuadEMfield* cpp_QuadEMfield = (QuadEMfield*)((pyORBIT_Object*) self)->cpp_obj; double czBy; //NO NEW OBJECT CREATED By0 PyArg_ParseTuple! NO NEED OF Py_DECREF() if(!PyArg_ParseTuple( args,"d:",&czBy)) error(" SetupPrint(Bz0 - parameter is needed"); - else + else cpp_QuadEMfield->czBy(czBy); - + Py_INCREF(Py_None); return Py_None; } - + static PyObject* QuadEMfield_cxBz(PyObject *self, PyObject *args){ - QuadEMfield* cpp_QuadEMfield = (QuadEMfield*)((pyORBIT_Object*) self)->cpp_obj; + QuadEMfield* cpp_QuadEMfield = (QuadEMfield*)((pyORBIT_Object*) self)->cpp_obj; double cxBz; //NO NEW OBJECT CREATED By0 PyArg_ParseTuple! NO NEED OF Py_DECREF() if(!PyArg_ParseTuple( args,"d:",&cxBz)) error(" SetupPrint(Bz0 - parameter is needed"); - else + else cpp_QuadEMfield->cxBz(cxBz); - + Py_INCREF(Py_None); return Py_None; } - + static PyObject* QuadEMfield_cyBz(PyObject *self, PyObject *args){ - QuadEMfield* cpp_QuadEMfield = (QuadEMfield*)((pyORBIT_Object*) self)->cpp_obj; + QuadEMfield* cpp_QuadEMfield = (QuadEMfield*)((pyORBIT_Object*) self)->cpp_obj; double cyBz; //NO NEW OBJECT CREATED By0 PyArg_ParseTuple! NO NEED OF Py_DECREF() if(!PyArg_ParseTuple( args,"d:",&cyBz)) error(" SetupPrint(Bz0 - parameter is needed"); - else + else cpp_QuadEMfield->cyBz(cyBz); - + Py_INCREF(Py_None); return Py_None; } - + static PyObject* QuadEMfield_czBz(PyObject *self, PyObject *args){ - QuadEMfield* cpp_QuadEMfield = (QuadEMfield*)((pyORBIT_Object*) self)->cpp_obj; + QuadEMfield* cpp_QuadEMfield = (QuadEMfield*)((pyORBIT_Object*) self)->cpp_obj; double czBz; //NO NEW OBJECT CREATED By0 PyArg_ParseTuple! NO NEED OF Py_DECREF() if(!PyArg_ParseTuple( args,"d:",&czBz)) error(" SetupPrint(Bz0 - parameter is needed"); - else + else cpp_QuadEMfield->czBz(czBz); - + Py_INCREF(Py_None); return Py_None; } - - - - - - - - - - + + + + + + + + + + //----------------------------------------------------- @@ -416,7 +414,7 @@ extern "C" { static void QuadEMfield_del(pyORBIT_Object* self){ //std::cerr<<"The LSFieldSource __del__ has been called!"<cpp_obj); - self->ob_type->tp_free((PyObject*)self); + Py_TYPE(self)->tp_free((PyObject*)self); } // defenition of the methods of the python PyBaseFieldSource wrapper class @@ -460,8 +458,7 @@ extern "C" { //new python PyBaseFieldSource wrapper type definition static PyTypeObject pyORBIT_QuadEMfield_Type = { - PyObject_HEAD_INIT(NULL) - 0, /*ob_size*/ + PyVarObject_HEAD_INIT(NULL, 0) "QuadEMfield", /*tp_name*/ sizeof(pyORBIT_Object), /*tp_basicsize*/ 0, /*tp_itemsize*/ @@ -499,7 +496,7 @@ extern "C" { (initproc) QuadEMfield_init, /* tp_init */ 0, /* tp_alloc */ QuadEMfield_new, /* tp_new */ - }; + }; //-------------------------------------------------- //Initialization function of the pyPyBaseFieldSource class diff --git a/README.md b/README.md index 32c6a83..4960a8b 100644 --- a/README.md +++ b/README.md @@ -1,34 +1,40 @@ [![Build Status](https://travis-ci.org/PyORBIT-Collaboration/laser-stripping.svg?branch=master)](https://travis-ci.org/PyORBIT-Collaboration/laser-stripping) -# Laser stripping extension for pyORBIT code +# Laser stripping extension for PyORBIT3 code Installation procedure requires building from source. -All installation steps happen in command line (terminal). +All installation steps happen in command line (terminal). - ## 1. Build pyORBIT as described [here](https://github.com/PyORBIT-Collaboration/py-orbit) + ## 1. Build PyORBIT3 as described [here](https://github.com/PyORBIT-Collaboration/PyORBIT3) ## 2. Clone the source code ```shell git clone https://github.com/PyORBIT-Collaboration/laser-stripping.git ``` -Put the cloned repository somewhere outside of pyORBIT directory. +Put the cloned repository somewhere outside of PyORBIT directory. So your source is now in *laser-stripping* directory. -## 3. Setup environment variables -If you built pyORBIT environment from source, use *customEnvironment.sh* instead. +## 3. Build & Install + +To compile with Meson + ```shell -source /setupEnvironment.sh +meson setup build --prefix=$CONDA_PREFIX # or your preferred path +meson compile -C build +meson install -C build ``` -## 4. Build the code -```shell -make clean -make -``` -This will put a dynamic library into *\/lib*. -If make failed, it usually means that some of the libraries aren't set up properly. +## 4. Running Examples -# Running Examples - ```shell cd examples -./START.sh 1GeV-particle-stripping.py 2 +python pyorbit3-1GeV-particle-stripping.py +# or, assuming PyORBIT was built with MPI enabled: +mpirun -np 2 python pyorbit3-1GeV-particle-stripping +``` +This will launch *pyorbit3-1GeV-particle-stripping.py* example on two MPI nodes. Other laser stripping related examples are availabale in [Examples](https://github.com/PyORBIT-Collaboration/examples/tree/master/ext/LaserStripping) repository. + +## 5. Uninstalling + +```shell +cd build +ninja uninstall +ninja clean ``` -This will launch *1GeV-particle-stripping.py* example on two MPI nodes. Other laser stripping related examples are availabale in [Examples](https://github.com/PyORBIT-Collaboration/examples/tree/master/ext/LaserStripping) repository. diff --git a/Record/PrintExtEffects.cc b/Record/PrintExtEffects.cc index 3b469e9..12d12f4 100644 --- a/Record/PrintExtEffects.cc +++ b/Record/PrintExtEffects.cc @@ -9,14 +9,14 @@ // DESCRIPTION // /////////////////////////////////////////////////////////////////////////// -#include "orbit_mpi.hh" +#include #include #include #include #include #include "PrintExtEffects.hh" -#include "RungeKuttaTracker.hh" +#include using namespace LaserStripping; using namespace OrbitUtils; @@ -27,48 +27,48 @@ PrintExtEffects::PrintExtEffects(std::string name,int i,std::string addr) setRankPrepare(1); setRankFinalize(1); setRankApply(1); - - + + ORBIT_MPI_Comm_size(MPI_COMM_WORLD, &size_MPI); ORBIT_MPI_Comm_rank(MPI_COMM_WORLD, &rank_MPI); - + num_print=i; addr_print=addr; eff_name = name; } -PrintExtEffects::~PrintExtEffects() +PrintExtEffects::~PrintExtEffects() { } -void PrintExtEffects::setupEffects(Bunch* bunch){ +void PrintExtEffects::setupEffects(Bunch* bunch){ setup_par = true; } - -void PrintExtEffects::applyEffectsForEach(Bunch* bunch, int i, - double* y_in_vct, double* y_out_vct, - double t, double t_step, + +void PrintExtEffects::applyEffectsForEach(Bunch* bunch, int i, + double* y_in_vct, double* y_out_vct, + double t, double t_step, OrbitUtils::BaseFieldSource* fieldSource, RungeKuttaTracker* tracker) { if (setup_par == true) { Num = tracker->getStepsNumber(); - + if(NumgetParticleAttributes(eff_name)->getAttSize();j++) + for (int j=0; jgetParticleAttributes(eff_name)->getAttSize();j++) { file<<"\t"<getParticleAttributes(eff_name)->attArr(i)[j]; } @@ -77,10 +77,10 @@ void PrintExtEffects::applyEffectsForEach(Bunch* bunch, int i, } } - - - - + + + + diff --git a/Record/PrintExtEffects.hh b/Record/PrintExtEffects.hh index 6a6d2e4..6592a6e 100644 --- a/Record/PrintExtEffects.hh +++ b/Record/PrintExtEffects.hh @@ -3,7 +3,7 @@ #include "Python.h" -#include "ExternalEffects.hh" +#include #include "BaseLaserFieldSource.hh" #include diff --git a/Record/RecordEvolution.cc b/Record/RecordEvolution.cc index 0deb064..a9f55f6 100644 --- a/Record/RecordEvolution.cc +++ b/Record/RecordEvolution.cc @@ -7,33 +7,33 @@ // 06/27/2008 // // DESCRIPTION -// The base class for C++ implementation of a external effects -// during the transport of particles through the external field. +// The base class for C++ implementation of a external effects +// during the transport of particles through the external field. // It should be sub-classed on Python level and implement // setupEffects(Bunch* bunch) // finalizeEffects(Bunch* bunch) -// applyEffects(Bunch* bunch, int index, -// double* y_in_vct, double* y_out_vct, -// double t, double t_step, +// applyEffects(Bunch* bunch, int index, +// double* y_in_vct, double* y_out_vct, +// double t, double t_step, // OrbitUtils::BaseFieldSource* fieldSource) // methods. // The results of these methods will be available from the c++ level. // This is an example of embedding Python in C++ Orbit level. // -// -// -// -// +// +// +// +// // /////////////////////////////////////////////////////////////////////////// -#include "orbit_mpi.hh" +#include #include #include #include #include #include "RecordEvolution.hh" -#include "RungeKuttaTracker.hh" +#include @@ -50,7 +50,7 @@ using namespace OrbitUtils; RecordEvolution::RecordEvolution(std::string effect,int ind_effect, int num) { - + setRankSetup(1); setRankPrepare(1); @@ -64,48 +64,48 @@ RecordEvolution::RecordEvolution(std::string effect,int ind_effect, int num) } -RecordEvolution::~RecordEvolution() +RecordEvolution::~RecordEvolution() { } void RecordEvolution::setupEffects(Bunch* bunch){ setup_par = true; - + if(bunch->hasParticleAttributes("Evolution")==1) bunch->removeParticleAttributes("Evolution"); - + std::map part_attr_dict; part_attr_dict["size"] = num_plot+5; bunch->addParticleAttributes("Evolution",part_attr_dict); - + Evol = bunch->getParticleAttributes("Evolution"); RecEff = bunch->getParticleAttributes(effect_name); - - for (int i=0; igetSize();i++) + + for (int i=0; igetSize();i++) Evol->attArr(i)[0] = RecEff->attArr(i)[index_effect]; } -void RecordEvolution::applyEffectsForEach(Bunch* bunch, int i, - double* y_in_vct, double* y_out_vct, - double t, double t_step, +void RecordEvolution::applyEffectsForEach(Bunch* bunch, int i, + double* y_in_vct, double* y_out_vct, + double t, double t_step, OrbitUtils::BaseFieldSource* fieldSource, RungeKuttaTracker* tracker) { - + if (setup_par == true) { - + x0 = y_in_vct[0];y0 = y_in_vct[1];z0 = y_in_vct[2]; Num = tracker->getStepsNumber(); - + if(NumgetSize();i++) { x0_ev(i) = x0; y0_ev(i) = y0; @@ -123,5 +123,5 @@ void RecordEvolution::applyEffectsForEach(Bunch* bunch, int i, } } - + } diff --git a/Record/RecordEvolution.hh b/Record/RecordEvolution.hh index 154c44f..8b8fd01 100644 --- a/Record/RecordEvolution.hh +++ b/Record/RecordEvolution.hh @@ -3,9 +3,9 @@ #include "Python.h" -#include "ExternalEffects.hh" +#include #include "BaseLaserFieldSource.hh" -#include "ParticleAttributes.hh" +#include #include using namespace TrackerRK4; diff --git a/Record/wrap_print_ext_effects.cc b/Record/wrap_print_ext_effects.cc index d570424..cfe4239 100644 --- a/Record/wrap_print_ext_effects.cc +++ b/Record/wrap_print_ext_effects.cc @@ -1,7 +1,6 @@ -#include "orbit_mpi.hh" -#include "pyORBIT_Object.hh" +#include +#include -#include "wrap_utils.hh" #include #include @@ -44,7 +43,7 @@ extern "C" { if(!PyArg_ParseTuple( args,"sis:",&eff_name,&max_print,&address)){ error("PrintExtEffects(eff_name,max_print,address) - params. are needed"); - } + } else { std::string name(eff_name); std::string addr(address); @@ -53,8 +52,8 @@ extern "C" { } return 0; } - - // name([name]) - sets or returns the name of the External Effeects class + + // name([name]) - sets or returns the name of the External Effeects class static PyObject* PrintExtEffects_name(PyObject *self, PyObject *args){ PrintExtEffects* cpp_PrintExtEffects = (PrintExtEffects*) ((pyORBIT_Object*) self)->cpp_obj; const char* name = NULL; @@ -66,8 +65,8 @@ extern "C" { cpp_PrintExtEffects->setName(name_str); } return Py_BuildValue("s",cpp_PrintExtEffects->getName().c_str()); - } - + } + //----------------------------------------------------- //destructor for python PyExternalEffects class (__del__ method). @@ -75,7 +74,7 @@ extern "C" { static void PrintExtEffects_del(pyORBIT_Object* self){ //std::cerr<<"The LasStripExternalEffects __del__ has been called!"<cpp_obj); - self->ob_type->tp_free((PyObject*)self); + Py_TYPE(self)->tp_free((PyObject*)self); } // defenition of the methods of the python PyExternalEffects wrapper class @@ -93,8 +92,7 @@ extern "C" { //new python PyExternalEffects wrapper type definition static PyTypeObject pyORBIT_PrintExtEffects_Type = { - PyObject_HEAD_INIT(NULL) - 0, /*ob_size*/ + PyVarObject_HEAD_INIT(NULL, 0) "PrintExtEffects", /*tp_name*/ sizeof(pyORBIT_Object), /*tp_basicsize*/ 0, /*tp_itemsize*/ @@ -132,10 +130,10 @@ extern "C" { (initproc) PrintExtEffects_init, /* tp_init */ 0, /* tp_alloc */ PrintExtEffects_new, /* tp_new */ - }; + }; + - //-------------------------------------------------- //Initialization function of the pyPyExternalEffects class //It will be called from Bunch wrapper initialization @@ -144,7 +142,7 @@ extern "C" { if (PyType_Ready(&pyORBIT_PrintExtEffects_Type) < 0) return; Py_INCREF(&pyORBIT_PrintExtEffects_Type); PyModule_AddObject(module, "PrintExtEffects", (PyObject *)&pyORBIT_PrintExtEffects_Type); - + } #ifdef __cplusplus diff --git a/Record/wrap_record_evolution.cc b/Record/wrap_record_evolution.cc index 6bb0ed3..1a924bc 100644 --- a/Record/wrap_record_evolution.cc +++ b/Record/wrap_record_evolution.cc @@ -1,7 +1,6 @@ -#include "orbit_mpi.hh" -#include "pyORBIT_Object.hh" +#include +#include -#include "wrap_utils.hh" #include @@ -43,23 +42,23 @@ extern "C" { int ind_effect; int num; PyObject* pyBaseLaserField=NULL; - + if(!PyArg_ParseTuple( args,"sii:",&effect, &ind_effect, &num)){ error("RecordEvolution(LaserField,delta_E,dipole_transition) - params. are needed"); - } + } else { BaseLaserFieldSource* lfs = (BaseLaserFieldSource*) ((pyORBIT_Object*) pyBaseLaserField)->cpp_obj; std::string ef(effect); self->cpp_obj = new RecordEvolution(ef,ind_effect, num); ((RecordEvolution*) self->cpp_obj)->setPyWrapper((PyObject*) self); } - + return 0; } - - - - // name([name]) - sets or returns the name of the External Effeects class + + + + // name([name]) - sets or returns the name of the External Effeects class static PyObject* RecordEvolution_name(PyObject *self, PyObject *args){ RecordEvolution* cpp_RecordEvolution = (RecordEvolution*) ((pyORBIT_Object*) self)->cpp_obj; const char* name = NULL; @@ -71,17 +70,17 @@ extern "C" { cpp_RecordEvolution->setName(name_str); } return Py_BuildValue("s",cpp_RecordEvolution->getName().c_str()); - } - - - + } + + + //----------------------------------------------------- //destructor for python PyExternalEffects class (__del__ method). //----------------------------------------------------- static void RecordEvolution_del(pyORBIT_Object* self){ //std::cerr<<"The LasStripExternalEffects __del__ has been called!"<cpp_obj); - self->ob_type->tp_free((PyObject*)self); + Py_TYPE(self)->tp_free((PyObject*)self); } // defenition of the methods of the python PyExternalEffects wrapper class @@ -99,8 +98,7 @@ extern "C" { //new python PyExternalEffects wrapper type definition static PyTypeObject pyORBIT_RecordEvolution_Type = { - PyObject_HEAD_INIT(NULL) - 0, /*ob_size*/ + PyVarObject_HEAD_INIT(NULL, 0) "RecordEvolution", /*tp_name*/ sizeof(pyORBIT_Object), /*tp_basicsize*/ 0, /*tp_itemsize*/ @@ -138,10 +136,10 @@ extern "C" { (initproc) RecordEvolution_init, /* tp_init */ 0, /* tp_alloc */ RecordEvolution_new, /* tp_new */ - }; + }; + - //-------------------------------------------------- //Initialization function of the pyPyExternalEffects class //It will be called from Bunch wrapper initialization @@ -150,7 +148,7 @@ extern "C" { if (PyType_Ready(&pyORBIT_RecordEvolution_Type) < 0) return; Py_INCREF(&pyORBIT_RecordEvolution_Type); PyModule_AddObject(module, "RecordEvolution", (PyObject *)&pyORBIT_RecordEvolution_Type); - + } #ifdef __cplusplus diff --git a/RegularGridFS/RegularGridFS.cc b/RegularGridFS/RegularGridFS.cc index 37df94c..6c27ba3 100644 --- a/RegularGridFS/RegularGridFS.cc +++ b/RegularGridFS/RegularGridFS.cc @@ -7,8 +7,8 @@ // 04/21/2003 // // DESCRIPTION -// The base class for Python implementation of a field source. -// It should be sub-classed on Python level and implements +// The base class for Python implementation of a field source. +// It should be sub-classed on Python level and implements // getElectricField(x,y,z,t) and getMagneticField (x,y,z,t) methods. // The results of these methods will be available from the c++ level. // This is an example of embedding Python in C++ Orbit level. @@ -17,8 +17,8 @@ /////////////////////////////////////////////////////////////////////////// #include "RegularGridFS.hh" -#include "orbit_mpi.hh" -#include "BufferStore.hh" +#include +#include #include #include @@ -30,90 +30,90 @@ using namespace OrbitUtils; RegularGridFS::RegularGridFS(std::string file_name,std::string type_of_field, double coord_mult, double field_mul) { orient = new FieldOrientation(); - + if (type_of_field[0]=='e') e_field=true; if (type_of_field[0]=='m') e_field=false; - + xyz_mult=coord_mult; field_mult = field_mul; - - + + int rank_MPI,size_MPI; ORBIT_MPI_Comm_size(MPI_COMM_WORLD, &size_MPI); ORBIT_MPI_Comm_rank(MPI_COMM_WORLD, &rank_MPI); - + std::ifstream file; double F,x,y,z; int fi; - + if(rank_MPI == 0) { file.open(file_name.c_str()); file>>x_min>>y_min>>z_min>>F>>F>>F; file.close(); - + file.open(file_name.c_str()); while(true) {file>>x>>y>>z>>F>>F>>F;if (x>x_min+fabs(x_min)*1.0e-14) {delta_x=x-x_min; break;}} file.close(); file.open(file_name.c_str()); while(true) {file>>x>>y>>z>>F>>F>>F;if (y>y_min+fabs(y_min)*1.0e-14) {delta_y=y-y_min; break;}} file.close(); - file.open(file_name.c_str()); while(true) {file>>x>>y>>z>>F>>F>>F;if (z>z_min+fabs(z_min)*1.0e-14) {delta_z=z-z_min; break;}} file.close(); - + file.open(file_name.c_str()); while(true) {file>>x>>y>>z>>F>>F>>F;if (z>z_min+fabs(z_min)*1.0e-14) {delta_z=z-z_min; break;}} file.close(); + file.open(file_name.c_str());fi=0; while(!file.eof()) {file>>x_max>>y_max>>z_max>>F>>F>>F;fi++;} file.close(); n_data=fi-1; - - + + x_min *=xyz_mult; y_min *=xyz_mult; z_min *=xyz_mult; x_max *=xyz_mult; y_max *=xyz_mult; z_max *=xyz_mult; delta_x *=xyz_mult; delta_y *=xyz_mult; delta_z *=xyz_mult; } - ORBIT_MPI_Bcast(&n_data, 1, MPI_INT,0,MPI_COMM_WORLD); - + ORBIT_MPI_Bcast(&n_data, 1, MPI_INT,0,MPI_COMM_WORLD); + ORBIT_MPI_Bcast(&x_min, 1, MPI_DOUBLE,0,MPI_COMM_WORLD); ORBIT_MPI_Bcast(&y_min, 1, MPI_DOUBLE,0,MPI_COMM_WORLD); ORBIT_MPI_Bcast(&z_min, 1, MPI_DOUBLE,0,MPI_COMM_WORLD); - + ORBIT_MPI_Bcast(&x_max, 1, MPI_DOUBLE,0,MPI_COMM_WORLD); ORBIT_MPI_Bcast(&y_max, 1, MPI_DOUBLE,0,MPI_COMM_WORLD); ORBIT_MPI_Bcast(&z_max, 1, MPI_DOUBLE,0,MPI_COMM_WORLD); - + ORBIT_MPI_Bcast(&delta_x, 1, MPI_DOUBLE,0,MPI_COMM_WORLD); ORBIT_MPI_Bcast(&delta_y, 1, MPI_DOUBLE,0,MPI_COMM_WORLD); ORBIT_MPI_Bcast(&delta_z, 1, MPI_DOUBLE,0,MPI_COMM_WORLD); - + nx=int((x_max-x_min)/delta_x+0.5)+1; ny=int((y_max-y_min)/delta_y+0.5)+1; nz=int((z_max-z_min)/delta_z+0.5)+1; - + nx_min = int(x_min/delta_x*(1+0.5/fabs(x_min/delta_x))); ny_min = int(y_min/delta_y*(1+0.5/fabs(y_min/delta_y))); nz_min = int(z_min/delta_z*(1+0.5/fabs(z_min/delta_z))); - - - B_x=new double**[nx+1]; - for (int i=0;igetFreeDoubleArr(buff_index,nz); - + //this loop reads magnetic fields if(rank_MPI == 0) { - file.open(file_name.c_str()); + file.open(file_name.c_str()); for(int i=0;isetUnusedDoubleArr(buff_index); + + + BufferStore::getBufferStore()->setUnusedDoubleArr(buff_index); } @@ -180,7 +180,7 @@ RegularGridFS::RegularGridFS(std::string file_name,std::string type_of_field, do RegularGridFS::~RegularGridFS() -{ +{ for (int i=0;isetCoefficients(x0, y0, z0,kx, ky, kz, mx, my, mz); - return; + orient->setCoefficients(x0, y0, z0,kx, ky, kz, mx, my, mz); + return; } double RegularGridFS::getNonOrientedBn(double*** B_n,double x, double y, double z, double t){ - + if ((x_min<= x&&x<= x_max)&&(y_min<= y&&y <= y_max)&&(z_min<= z&&z <= z_max)) { - - + + double x_t = x/delta_x; double y_t = y/delta_y; double z_t = z/delta_z; - + int nx = int(x_t); int ny = int(y_t); int nz = int(z_t); - -/* + +/* for (int i=0; i<2; i++) for (int j=0; j<2; j++) for (int k=0; k<2; k++) - + sum += pow(-1,i+j+k)*B_n[nx-nx_min+i][ny-ny_min+j][nz-nz_min+k]*(x_t-(nx+1-i))*(y_t-(ny+1-j))*(z_t-(nz+1-k)); -*/ return +*/ return -B_n[nx-nx_min][ny-ny_min][nz-nz_min]*(x_t-(nx+1))*(y_t-(ny+1))*(z_t-(nz+1)) +B_n[nx-nx_min][ny-ny_min][nz-nz_min+1]*(x_t-(nx+1))*(y_t-(ny+1))*(z_t-nz) @@ -232,64 +232,64 @@ double RegularGridFS::getNonOrientedBn(double*** B_n,double x, double y, double +B_n[nx-nx_min+1][ny-ny_min+1][nz-nz_min+1]*(x_t-nx)*(y_t-ny)*(z_t-nz); } - + else return 0; - - - - + + + + } -void RegularGridFS::getConstructedNonOrientedField(double& H_x,double& H_y,double& H_z, - double*** B_x,double*** B_y,double*** B_z, +void RegularGridFS::getConstructedNonOrientedField(double& H_x,double& H_y,double& H_z, + double*** B_x,double*** B_y,double*** B_z, double x,double y,double z,double& t){ - + if (z>=0) { - + H_x=getNonOrientedBn(B_x,x,y,z,t); H_y=getNonOrientedBn(B_y,x,y,z,t); H_z=getNonOrientedBn(B_z,x,y,z,t); - + } - + if (z<0) { - + H_x = getNonOrientedBn(B_x,x,y,-z,t); H_y = getNonOrientedBn(B_y,x,y,-z,t); H_z =-getNonOrientedBn(B_z,x,y,-z,t); - + } - - + + return; - + } - -void RegularGridFS::getElectricMagneticField(double x, double y, double z, double t, + +void RegularGridFS::getElectricMagneticField(double x, double y, double z, double t, double& E_x, double& E_y, double& E_z, double& H_x, double& H_y, double& H_z) { - + double F_x, F_y, F_z; orient->OrientCoordinates(x,y,z); - getConstructedNonOrientedField(F_x, F_y, F_z, B_x, B_y, B_z, x, y, z, t); + getConstructedNonOrientedField(F_x, F_y, F_z, B_x, B_y, B_z, x, y, z, t); orient->OrientVector(F_x,F_y,F_z); - - + + if (e_field) { H_x = 0.0; H_y = 0.0; H_z = 0.0; E_x = F_x; E_y = F_y; E_z = F_z;} if (!e_field) { E_x = 0.0; E_y = 0.0; E_z = 0.0; H_x = F_x; H_y = F_y; H_z = F_z;} - -// std::cout<< e_field< #include "FieldOrientation.hh" #include #include diff --git a/RegularGridFS/wrap_regular_grid_fs.cc b/RegularGridFS/wrap_regular_grid_fs.cc index fc77271..e0b2004 100644 --- a/RegularGridFS/wrap_regular_grid_fs.cc +++ b/RegularGridFS/wrap_regular_grid_fs.cc @@ -1,7 +1,6 @@ -#include "orbit_mpi.hh" -#include "pyORBIT_Object.hh" +#include +#include -#include "wrap_utils.hh" #include "wrap_regular_grid_fs.hh" #include @@ -10,7 +9,6 @@ #include "RegularGridFS.hh" using namespace OrbitUtils; -using namespace wrap_orbit_utils; namespace wrap_regular_grid_fs{ @@ -36,43 +34,43 @@ extern "C" { //initializator for python CppBaseFieldSource class //this is implementation of the __init__ method - - + + static int RegularGridFS_init(pyORBIT_Object *self, PyObject *args, PyObject *kwds){ - - + + const char* file_name; const char* type_of_field; double field_mult; double coord_mult; - + if(!PyArg_ParseTuple( args,"ssdd:",&file_name, &type_of_field, &coord_mult, &field_mult)){ error("HydrogenStarkParam(file_name , type_of_field, coodr_mult, field_mult) - params. are needed"); - } + } else { - + std::string name_str_name(file_name); std::string name_str_field(type_of_field); self->cpp_obj = new RegularGridFS(name_str_name, name_str_field, coord_mult, field_mult); ((RegularGridFS*) self->cpp_obj)->setPyWrapper((PyObject*) self); } - - + + return 0; } - - - - + + + + static PyObject* RegularGridFS_setFieldOrientation(PyObject *self, PyObject *args){ RegularGridFS* RegGrid = (RegularGridFS*)((pyORBIT_Object*) self)->cpp_obj; - + double x0; double y0; double z0; @@ -87,27 +85,27 @@ extern "C" { if(!PyArg_ParseTuple( args,"ddddddddd:",&x0, &y0, &z0, &kx, &ky, &kz, &mx, &my, &mz)) {error("LaserExternalEfects - setLaserHalfAngle(x0, y0, z0, kx, ky, kz, mx, my, mz) - params. are needed");flag=false;} - + if(fabs(kx*mx+ky*my+kz*mz)/sqrt(kx*kx+ky*ky+kz*kz)/sqrt(mx*mx+my*my+mz*mz)>1.e-15) {error("Please be shure that kx*mx+ky*my+kz*mz==0");flag=false;} - - if(flag) + + if(flag) RegGrid->setFieldOrientation(x0,y0, z0, kx, ky, kz, mx, my, mz); - - + + Py_INCREF(Py_None); - return Py_None; - } - + return Py_None; + } + + + - - static PyObject* RegularGridFS_getFields(PyObject *self, PyObject *args){ RegularGridFS* cpp_fields = (RegularGridFS*)((pyORBIT_Object*) self)->cpp_obj; - - + + int nVars = PyTuple_Size(args); - + double x; double y; double z; @@ -118,19 +116,19 @@ extern "C" { double Bx; double By; double Bz; - + //NO NEW OBJECT CREATED BY PyArg_ParseTuple! NO NEED OF Py_DECREF() if(!PyArg_ParseTuple( args,"dddd:",&x,&y,&z,&t)) error(" getFields(x,y,z,t) - parameters are needed"); else cpp_fields->getElectricMagneticField(x,y,z,t,Ex,Ey,Ez,Bx,By,Bz); - + if (Ex==0.0&&Ey==0.0&&Ez==0.0) return Py_BuildValue("ddd",Bx,By,Bz); if (Bx==0.0&&By==0.0&&Bz==0.0) return Py_BuildValue("ddd",Ex,Ey,Ez); } - - + + //----------------------------------------------------- @@ -139,7 +137,7 @@ extern "C" { static void RegularGridFS_del(pyORBIT_Object* self){ //std::cerr<<"The RegularGridFS __del__ has been called!"<cpp_obj); - self->ob_type->tp_free((PyObject*)self); + Py_TYPE(self)->tp_free((PyObject*)self); } // defenition of the methods of the python PyBaseFieldSource wrapper class @@ -165,8 +163,7 @@ extern "C" { //new python PyBaseFieldSource wrapper type definition static PyTypeObject pyORBIT_RegularGridFS_Type = { - PyObject_HEAD_INIT(NULL) - 0, /*ob_size*/ + PyVarObject_HEAD_INIT(NULL, 0) "RegularGridFS", /*tp_name*/ sizeof(pyORBIT_Object), /*tp_basicsize*/ 0, /*tp_itemsize*/ @@ -204,7 +201,7 @@ extern "C" { (initproc) RegularGridFS_init, /* tp_init */ 0, /* tp_alloc */ RegularGridFS_new, /* tp_new */ - }; + }; //-------------------------------------------------- //Initialization function of the pyPyBaseFieldSource class diff --git a/SNSstrippingMagnet/SNSstrippingMagnet.cc b/SNSstrippingMagnet/SNSstrippingMagnet.cc index 8baa603..987879b 100644 --- a/SNSstrippingMagnet/SNSstrippingMagnet.cc +++ b/SNSstrippingMagnet/SNSstrippingMagnet.cc @@ -7,8 +7,8 @@ // 04/21/2003 // // DESCRIPTION -// The base class for Python implementation of a field source. -// It should be sub-classed on Python level and implements +// The base class for Python implementation of a field source. +// It should be sub-classed on Python level and implements // getElectricField(x,y,z,t) and getMagneticField (x,y,z,t) methods. // The results of these methods will be available from the c++ level. // This is an example of embedding Python in C++ Orbit level. @@ -17,7 +17,7 @@ /////////////////////////////////////////////////////////////////////////// //#include -#include "orbit_mpi.hh" +#include //#include ////#include //#include @@ -31,27 +31,27 @@ using namespace OrbitUtils; SNSstrippingMagnet::SNSstrippingMagnet(std::string field_data) -{ +{ pos = 0; len = 0; - std::ifstream file; + std::ifstream file; file.open(field_data.c_str()); - std::string line; + std::string line; while (std::getline(file,line)) { len++;} file.close(); - + z = new double [len]; Bx = new double[len]; - - file.open(field_data.c_str()); - for(int i = 0; i < len; i++){file>>z[i]>>Bx[i];} - file.close(); - + + file.open(field_data.c_str()); + for(int i = 0; i < len; i++){file>>z[i]>>Bx[i];} + file.close(); + A = 1; l = z[len- 1] - z[0]; - dz = l/(len - 1); + dz = l/(len - 1); } @@ -59,7 +59,7 @@ SNSstrippingMagnet::SNSstrippingMagnet(std::string field_data) SNSstrippingMagnet::~SNSstrippingMagnet() -{ +{ delete [] z; delete [] Bx; @@ -69,11 +69,11 @@ SNSstrippingMagnet::~SNSstrippingMagnet() double SNSstrippingMagnet::getBx(double zi){ - - + + if (zi <= z[0] || zi >= z[len - 1]) {return 0; } else{ - int i = (int)((zi - z[0])/dz); + int i = (int)((zi - z[0])/dz); return Bx[i] + (Bx[i+1] - Bx[i])*(zi - z[i])/dz; } } @@ -83,12 +83,12 @@ double SNSstrippingMagnet::getBx(double zi){ void SNSstrippingMagnet::setA(double _A){A = _A;} void SNSstrippingMagnet::setPosition(double _pos){ pos = _pos;} -void SNSstrippingMagnet::getElectricMagneticField(double x, double y, double z, double t, +void SNSstrippingMagnet::getElectricMagneticField(double x, double y, double z, double t, double& E_x, double& E_y, double& E_z, double& B_x, double& B_y, double& B_z) -{ - +{ + B_x = getBx(z - pos); B_y = 0; B_z = 0; @@ -105,4 +105,4 @@ return; - \ No newline at end of file + diff --git a/SNSstrippingMagnet/SNSstrippingMagnet.hh b/SNSstrippingMagnet/SNSstrippingMagnet.hh index bdcf451..2160606 100644 --- a/SNSstrippingMagnet/SNSstrippingMagnet.hh +++ b/SNSstrippingMagnet/SNSstrippingMagnet.hh @@ -24,7 +24,7 @@ #define SNSSTRIPPINGMAGNET_HH_ #include "Python.h" -#include "BaseFieldSource.hh" +#include #include diff --git a/SNSstrippingMagnet/wrap_snsstrippingmagnet.cc b/SNSstrippingMagnet/wrap_snsstrippingmagnet.cc index 50ca54e..c4e8632 100644 --- a/SNSstrippingMagnet/wrap_snsstrippingmagnet.cc +++ b/SNSstrippingMagnet/wrap_snsstrippingmagnet.cc @@ -1,7 +1,6 @@ -#include "orbit_mpi.hh" -#include "pyORBIT_Object.hh" +#include +#include -#include "wrap_utils.hh" #include "wrap_snsstrippingmagnet.hh" #include @@ -10,7 +9,6 @@ #include "SNSstrippingMagnet.hh" using namespace OrbitUtils; -using namespace wrap_orbit_utils; namespace wrap_snsstrippingmagnet{ @@ -36,34 +34,34 @@ extern "C" { //initializator for python CppBaseFieldSource class //this is implementation of the __init__ method - - + + static int SNSstrippingMagnet_init(pyORBIT_Object *self, PyObject *args, PyObject *kwds){ - - - const char* file_name; + + + const char* file_name; if(!PyArg_ParseTuple( args,"s:",&file_name)){ error("SNSstrippingMagnet(file_name ) - params. are needed"); - } + } else { - + std::string name_str_name(file_name); self->cpp_obj = new SNSstrippingMagnet(name_str_name); ((SNSstrippingMagnet*) self->cpp_obj)->setPyWrapper((PyObject*) self); } - - + + return 0; } - - + + static PyObject* SNSstrippingMagnet_setA(PyObject *self, PyObject *args){ SNSstrippingMagnet* cpp_field = (SNSstrippingMagnet*)((pyORBIT_Object*) self)->cpp_obj; @@ -73,15 +71,15 @@ extern "C" { error(" setA(z) - parameters are needed"); else cpp_field->setA(A); - + Py_INCREF(Py_None); - return Py_None; - - } - + return Py_None; + + } + + + - - static PyObject* SNSstrippingMagnet_setPosition(PyObject *self, PyObject *args){ SNSstrippingMagnet* cpp_field = (SNSstrippingMagnet*)((pyORBIT_Object*) self)->cpp_obj; double pos; @@ -91,11 +89,11 @@ extern "C" { else cpp_field->setPosition(pos); Py_INCREF(Py_None); - return Py_None; - - } - - // name([name]) - sets or returns the name of the External Effeects class + return Py_None; + + } + + // name([name]) - sets or returns the name of the External Effeects class static PyObject* SNSstrippingMagnet_getBx(PyObject *self, PyObject *args){ SNSstrippingMagnet* cpp_field = (SNSstrippingMagnet*)((pyORBIT_Object*) self)->cpp_obj; @@ -107,7 +105,7 @@ extern "C" { return Py_BuildValue("d",cpp_field->getBx(z)); } - + //----------------------------------------------------- @@ -116,7 +114,7 @@ extern "C" { static void SNSstrippingMagnet_del(pyORBIT_Object* self){ //std::cerr<<"The SNSstrippingMagnet __del__ has been called!"<cpp_obj); - self->ob_type->tp_free((PyObject*)self); + Py_TYPE(self)->tp_free((PyObject*)self); } // defenition of the methods of the python PyBaseFieldSource wrapper class @@ -140,8 +138,7 @@ extern "C" { //new python PyBaseFieldSource wrapper type definition static PyTypeObject pyORBIT_SNSstrippingMagnet_Type = { - PyObject_HEAD_INIT(NULL) - 0, /*ob_size*/ + PyVarObject_HEAD_INIT(NULL, 0) "SNSstrippingMagnet", /*tp_name*/ sizeof(pyORBIT_Object), /*tp_basicsize*/ 0, /*tp_itemsize*/ @@ -179,7 +176,7 @@ extern "C" { (initproc) SNSstrippingMagnet_init, /* tp_init */ 0, /* tp_alloc */ SNSstrippingMagnet_new, /* tp_new */ - }; + }; //-------------------------------------------------- //Initialization function of the pyPyBaseFieldSource class diff --git a/SchrodingerEquation/SchrodingerEquation.cc b/SchrodingerEquation/SchrodingerEquation.cc index 23961fd..7838fe9 100644 --- a/SchrodingerEquation/SchrodingerEquation.cc +++ b/SchrodingerEquation/SchrodingerEquation.cc @@ -7,26 +7,26 @@ // 06/27/2008 // // DESCRIPTION -// The base class for C++ implementation of a external effects -// during the transport of particles through the external field. +// The base class for C++ implementation of a external effects +// during the transport of particles through the external field. // It should be sub-classed on Python level and implement // setupEffects(Bunch* bunch) // finalizeEffects(Bunch* bunch) -// applyEffects(Bunch* bunch, int index, -// double* y_in_vct, double* y_out_vct, -// double t, double t_step, +// applyEffects(Bunch* bunch, int index, +// double* y_in_vct, double* y_out_vct, +// double t, double t_step, // OrbitUtils::BaseFieldSource* fieldSource) // methods. // The results of these methods will be available from the c++ level. // This is an example of embedding Python in C++ Orbit level. // -// -// -// -// +// +// +// +// // /////////////////////////////////////////////////////////////////////////// -#include "orbit_mpi.hh" +#include #include #include #include @@ -35,9 +35,9 @@ #include "SchrodingerEquation.hh" -#include "RungeKuttaTracker.hh" -#include "OrbitConst.hh" -#include "LorentzTransformationEM.hh" +#include +#include +//#include "LorentzTransformationEM.hh" #include "FieldRotation.hh" @@ -57,21 +57,21 @@ using namespace OrbitUtils; SchrodingerEquation::SchrodingerEquation(BaseLaserFieldSource* BaseLaserField, Stark* Stark,double par_res) { setName("unnamed"); - + StarkEffect=Stark; LaserField=BaseLaserField; Parameter_resonance=par_res; levels = StarkEffect->getStates()*(1+StarkEffect->getStates())*(1+2*StarkEffect->getStates())/6; - - - + + + //allocating memory for koefficients of 4-th order Runge-Kutta method and other koeeficients of the master equation k_RungeKutt=new tcomplex*[levels+1]; for (int i=0;igetPyWrapper() == NULL){ delete LaserField; - } else { + } else { Py_XDECREF(LaserField->getPyWrapper()); } - - + + if(StarkEffect->getPyWrapper() == NULL){ - delete StarkEffect; + delete StarkEffect; } else { Py_XDECREF(StarkEffect->getPyWrapper()); } @@ -139,88 +139,88 @@ SchrodingerEquation::~SchrodingerEquation() void SchrodingerEquation::CalcPopulations(int i, Bunch* bunch) { - + PopAttr->attArr(i)[0] = 1; for (int j=1; jattArr(i)[j] = Re(i,j)*Re(i,j)+Im(i,j)*Im(i,j); PopAttr->attArr(i)[0] -= PopAttr->attArr(i)[j]; } - + } -void SchrodingerEquation::setupEffects(Bunch* bunch){ - - - - +void SchrodingerEquation::setupEffects(Bunch* bunch){ + + + + int_E=new tcomplex**[levels+1]; for (int i=0;igetSize()]; for(int i=0; igetSize();i++) for(int n=1; nhasParticleAttributes("Amplitudes")==0) { std::map part_attr_dict; part_attr_dict["size"] = 2*levels+1; bunch->addParticleAttributes("Amplitudes",part_attr_dict); - + for (int i=0; igetSize();i++) bunch->getParticleAttributes("Amplitudes")->attValue(i,1) = 1; } - + if(bunch->hasParticleAttributes("Populations")==0) { std::map part_attr_dict; part_attr_dict["size"] = levels+1; bunch->addParticleAttributes("Populations",part_attr_dict); } - - + + AmplAttr = bunch->getParticleAttributes("Amplitudes"); PopAttr = bunch->getParticleAttributes("Populations"); - - - + + + nx=new double[bunch->getSize()]; ny=new double[bunch->getSize()]; nz=new double[bunch->getSize()]; - + install_field_dir=new bool [bunch->getSize()]; - + for (int i=0; igetSize();i++) install_field_dir[i] = true; - - - + + + for (int i=0; igetSize();i++) CalcPopulations(i, bunch); } - + //-0.00226674 0.000171211 -0.326731 0.000842283-0.000358742 1.69601 - + void SchrodingerEquation::finalizeEffects(Bunch* bunch){ - + for (int i=0;igetElectricMagneticField(x0,y0,z0,t,Ex_stat,Ey_stat,Ez_stat,Bx_stat,By_stat,Bz_stat); - LorentzTransformationEM::transform(bunch->getMass(),px0,py0,pz0,Ex_stat,Ey_stat,Ez_stat,Bx_stat,By_stat,Bz_stat); - + + fieldSource->getElectricMagneticField(x0,y0,z0,t,Ex_stat,Ey_stat,Ez_stat,Bx_stat,By_stat,Bz_stat); + //LorentzTransformationEM::transform(bunch->getMass(),px0,py0,pz0,Ex_stat,Ey_stat,Ez_stat,Bx_stat,By_stat,Bz_stat); + if (install_field_dir[i]) { - nx[i] = Ex_stat; - ny[i] = Ey_stat; + nx[i] = Ex_stat; + ny[i] = Ey_stat; nz[i] = Ez_stat; install_field_dir[i]=false; - } - - + } + + Ez_stat=sqrt(Ex_stat*Ex_stat+Ey_stat*Ey_stat+Ez_stat*Ez_stat); - Ex_stat=0; + Ex_stat=0; Ey_stat=0; - + if(LaserField->region(x,y,z)) { - + for (int j=0; j<3;j++) { - + LaserField->getLaserEMField(x0+j*(x-x0)/2,y0+j*(y-y0)/2,z0+j*(z-z0)/2, t+j*t_step/2,Ex_las[j],Ey_las[j],Ez_las[j],Bx_las[j],By_las[j],Bz_las[j]); - - - LorentzTransformationEM::complex_transform(bunch->getMass(), - px0,py0,pz0, - Ex_las[j],Ey_las[j],Ez_las[j], - Bx_las[j],By_las[j],Bz_las[j]); + + + //LorentzTransformationEM::complex_transform(bunch->getMass(), + // px0,py0,pz0, + // Ex_las[j],Ey_las[j],Ez_las[j], + // Bx_las[j],By_las[j],Bz_las[j]); FieldRotation::RotateElectricFieldsV(nx[i],ny[i],nz[i],Ex_las[j],Ey_las[j],Ez_las[j]); - + } } - else - for(int j=0;j<3;j++) { + else + for(int j=0;j<3;j++) { Ex_las[j] = 0.; Ey_las[j] = 0.; Ez_las[j] = 0.; } - - + + } @@ -327,8 +327,8 @@ void SchrodingerEquation::GetParticleFrameFields(int i,double t, double t_step, void SchrodingerEquation::GetParticleFrameParameters(int i, double t,double t_step, Bunch* bunch) { - - + + double ta=2.418884326505e-17; //atomic unit of time double Ea=5.14220642e011; //Atomic unit of electric field double m=bunch->getMass(); @@ -342,21 +342,21 @@ part_t_step=t_step*coeff; //time step in frame of particle (in atomic units) -for(int j=0;j<3;j++) -{ - Ex_las[j]/=Ea; - Ey_las[j]/=Ea; +for(int j=0;j<3;j++) +{ + Ex_las[j]/=Ea; + Ey_las[j]/=Ea; Ez_las[j]/=Ea; } -Ex_stat/=Ea; -Ey_stat/=Ea; +Ex_stat/=Ea; +Ey_stat/=Ea; Ez_stat/=Ea; omega_part=gamma*ta*LaserField->getFrequencyOmega(m,x0,y0,z0,px0,py0,pz0,t); // frequensy of laser in particle frame (in atomic units) - + } @@ -368,52 +368,52 @@ omega_part=gamma*ta*LaserField->getFrequencyOmega(m,x0,y0,z0,px0,py0,pz0,t); // void SchrodingerEquation::AmplSolver4step(int i, Bunch* bunch) { - + tcomplex z,z1,z2,JdEdt, mu_x,mu_y,mu_z; double dt,dE; - - - StarkEffect->SetE(Ez_stat); - - - + StarkEffect->SetE(Ez_stat); + + + + + for(int n=2; nfield_thresh[n]>Ez_stat && StarkEffect->field_thresh[m]>Ez_stat) { - + StarkEffect->getTransition(n,m,mu_x,mu_y,mu_z); dE = fabs(StarkEffect->En[n]-StarkEffect->En[m]); JdEdt = tcomplex(0.,dE*part_t_step); - + mu_Elas[n][m][1]=(mu_x*conj(Ex_las[1])+mu_y*conj(Ey_las[1])+mu_z*conj(Ez_las[1]))*(J/2.); cond[n][m]=fabs(dE - omega_part)field_thresh[m]>Ez_stat) { - + if (j==4) dt=part_t_step; else dt=part_t_step/2.; k_RungeKutt[m][j]=tcomplex(0); - for(int k=1;kfield_thresh[k]>Ez_stat) k_RungeKutt[m][j]-=exp_mu_El[m][k][j]*(a(i,k)+k_RungeKutt[k][j-1]*dt); - for(int k=m+1;kfield_thresh[k]>Ez_stat) k_RungeKutt[m][j]+=conj(exp_mu_El[k][m][j])*(a(i,k)+k_RungeKutt[k][j-1]*dt); - - + for(int k=1;kfield_thresh[k]>Ez_stat) k_RungeKutt[m][j]-=exp_mu_El[m][k][j]*(a(i,k)+k_RungeKutt[k][j-1]*dt); + for(int k=m+1;kfield_thresh[k]>Ez_stat) k_RungeKutt[m][j]+=conj(exp_mu_El[k][m][j])*(a(i,k)+k_RungeKutt[k][j-1]*dt); + + k_RungeKutt[m][j]-=StarkEffect->Gamman[m]*(a(i,m)+k_RungeKutt[m][j-1]*dt)/2.; - + } - - - - - + + + + + for(int m=1;mfield_thresh[m]>Ez_stat) { - z=part_t_step*(k_RungeKutt[m][1]+2.*k_RungeKutt[m][2]+2.*k_RungeKutt[m][3]+k_RungeKutt[m][4])/6.; + z=part_t_step*(k_RungeKutt[m][1]+2.*k_RungeKutt[m][2]+2.*k_RungeKutt[m][3]+k_RungeKutt[m][4])/6.; Re(i,m)+=z.real(); Im(i,m)+=z.imag(); } @@ -453,9 +453,9 @@ void SchrodingerEquation::AmplSolver4step(int i, Bunch* bunch) { Re(i,m) = 0; Im(i,m) = 0; } - + } - + } diff --git a/SchrodingerEquation/SchrodingerEquation.hh b/SchrodingerEquation/SchrodingerEquation.hh index 9c533b9..aade679 100644 --- a/SchrodingerEquation/SchrodingerEquation.hh +++ b/SchrodingerEquation/SchrodingerEquation.hh @@ -5,10 +5,10 @@ #include "Python.h" -#include "ExternalEffects.hh" +#include #include "BaseLaserFieldSource.hh" #include "Stark.hh" -#include "ParticleAttributes.hh" +#include using namespace TrackerRK4; diff --git a/SchrodingerEquation/wrap_schrodinger_equation.cc b/SchrodingerEquation/wrap_schrodinger_equation.cc index dcd37a9..0fc4d23 100644 --- a/SchrodingerEquation/wrap_schrodinger_equation.cc +++ b/SchrodingerEquation/wrap_schrodinger_equation.cc @@ -1,7 +1,6 @@ -#include "orbit_mpi.hh" -#include "pyORBIT_Object.hh" +#include +#include -#include "wrap_utils.hh" #include "wrap_schrodinger_equation.hh" #include @@ -39,7 +38,7 @@ extern "C" { //initializator for python PyExternalEffects class //this is implementation of the __init__ method static int SchrodingerEquation_init(pyORBIT_Object *self, PyObject *args, PyObject *kwds){ - + double par_res=0; @@ -49,21 +48,21 @@ extern "C" { if(!PyArg_ParseTuple( args,"OOd:",&pyBaseLaserField,&pyStarkEffect,&par_res)){ error("SchrodingerEquation(LaserField,StarkEffect,par_res) - params. are needed"); - } + } else { BaseLaserFieldSource* lfs = (BaseLaserFieldSource*) ((pyORBIT_Object*) pyBaseLaserField)->cpp_obj; Stark* StarkEffect = (Stark*) ((pyORBIT_Object*) pyStarkEffect)->cpp_obj; self->cpp_obj = new SchrodingerEquation(lfs, StarkEffect, par_res); ((SchrodingerEquation*) self->cpp_obj)->setPyWrapper((PyObject*) self); } - + return 0; } - - - - // name([name]) - sets or returns the name of the External Effeects class + + + + // name([name]) - sets or returns the name of the External Effeects class static PyObject* SchrodingerEquation_name(PyObject *self, PyObject *args){ SchrodingerEquation* cpp_SchrodingerEquation = (SchrodingerEquation*) ((pyORBIT_Object*) self)->cpp_obj; const char* name = NULL; @@ -75,20 +74,20 @@ extern "C" { cpp_SchrodingerEquation->setName(name_str); } return Py_BuildValue("s",cpp_SchrodingerEquation->getName().c_str()); - } - - - - - - + } + + + + + + //----------------------------------------------------- //destructor for python PyExternalEffects class (__del__ method). //----------------------------------------------------- static void SchrodingerEquation_del(pyORBIT_Object* self){ //std::cerr<<"The DensityMatrix __del__ has been called!"<cpp_obj); - self->ob_type->tp_free((PyObject*)self); + Py_TYPE(self)->tp_free((PyObject*)self); } // defenition of the methods of the python PyExternalEffects wrapper class @@ -108,8 +107,7 @@ extern "C" { //new python PyExternalEffects wrapper type definition static PyTypeObject pyORBIT_SchrodingerEquation_Type = { - PyObject_HEAD_INIT(NULL) - 0, /*ob_size*/ + PyVarObject_HEAD_INIT(NULL, 0) "SchrodingerEquation", /*tp_name*/ sizeof(pyORBIT_Object), /*tp_basicsize*/ 0, /*tp_itemsize*/ @@ -147,10 +145,10 @@ extern "C" { (initproc) SchrodingerEquation_init, /* tp_init */ 0, /* tp_alloc */ SchrodingerEquation_new, /* tp_new */ - }; + }; + - //-------------------------------------------------- //Initialization function of the pyPyExternalEffects class //It will be called from Bunch wrapper initialization @@ -159,7 +157,7 @@ extern "C" { if (PyType_Ready(&pyORBIT_SchrodingerEquation_Type) < 0) return; Py_INCREF(&pyORBIT_SchrodingerEquation_Type); PyModule_AddObject(module, "SchrodingerEquation", (PyObject *)&pyORBIT_SchrodingerEquation_Type); - + } #ifdef __cplusplus diff --git a/Stark/Stark.cc b/Stark/Stark.cc index 138046d..864f92e 100644 --- a/Stark/Stark.cc +++ b/Stark/Stark.cc @@ -7,19 +7,19 @@ // 04/21/2003 // // DESCRIPTION -// The base class for Python implementation of a field source. -// It should be sub-classed on Python level and implements +// The base class for Python implementation of a field source. +// It should be sub-classed on Python level and implements // getElectricField(x,y,z,t) and getMagneticField (x,y,z,t) methods. // The results of these methods will be available from the c++ level. // This is an example of embedding Python in C++ Orbit level. // // /////////////////////////////////////////////////////////////////////////// -#include "orbit_mpi.hh" -#include "BufferStore.hh" +#include +#include #include "Stark.hh" -#include "LorentzTransformationEM.hh" +//#include "LorentzTransformationEM.hh" #include #include #include @@ -37,12 +37,12 @@ inline int convert3to1level(int n,int n1, int m){ Stark::Stark(std::string addressEG,int states) -{ +{ int rank_MPI,size_MPI; ORBIT_MPI_Comm_size(MPI_COMM_WORLD, &size_MPI); ORBIT_MPI_Comm_rank(MPI_COMM_WORLD, &rank_MPI); - + std::ifstream file; double F; char G_ch[200]; @@ -51,185 +51,185 @@ Stark::Stark(std::string addressEG,int states) char nameEG[MAX_LENGTH_ADDRESS]; std::string dump,dump2; - + st=states; const_relax = 4*alpha*alpha*alpha/3; levels=st*(1+st)*(1+2*st)/6; - + field_thresh = new double[levels+1]; n_data = new int[levels+1]; En = new double[levels+1]; - Gamman = new double[levels+1]; - + Gamman = new double[levels+1]; + + - if(rank_MPI == 0) { n_data[0] = 1; - - snprintf(nameEG,MAX_LENGTH_ADDRESS,"%s000.txt",addressEG.c_str()); + + snprintf(nameEG,MAX_LENGTH_ADDRESS,"%s000.txt",addressEG.c_str()); file.open(nameEG); file>>F>>F>>F>>delta_F; file.clear();file.close(); - + snprintf(nameEG,MAX_LENGTH_ADDRESS,"%stransitions.txt",addressEG.c_str()); file.open(nameEG); file>>ff>>ff>>ff>>dump>>ff>>ff>>ff; fi=0; while(dump2!=dump) {file>>dump2;fi++;} file.clear();file.close(); order_trans = (fi-4)/3; - - + + for(int n=1;n>field_thresh[k]>>F>>G_ch; n_data[k]++;} file.clear();file.close(); - + k=convert3to1level(n,n1,m); + snprintf(nameEG,MAX_LENGTH_ADDRESS,"%s%i%i%i.txt",addressEG.c_str(),n1,n-n1-abs(m)-1,abs(m)); + n_data[k]=-1; file.open(nameEG); while(!file.eof()) {file>>field_thresh[k]>>F>>G_ch; n_data[k]++;} file.clear();file.close(); + } } } - - } + + } + - - ORBIT_MPI_Bcast(&order_trans, 1, MPI_INT,0,MPI_COMM_WORLD); + ORBIT_MPI_Bcast(&order_trans, 1, MPI_INT,0,MPI_COMM_WORLD); ORBIT_MPI_Bcast(&delta_F, 1, MPI_DOUBLE,0,MPI_COMM_WORLD); for (int i=0;i>F>>energy[k][fi]>>G_ch; gamma_autoionization[k][fi] = atof(G_ch);} file.clear();file.close(); } - - + + for(int fi=0;fi>nn1>>nn2>>mm>>dump>>nn1s>>nn2s>>mms; - + if(nn1==n1&&nn2==n-n1-fabs(m)-1&&mm==m&&nn1s==n1s&&nn2s==ns-n1s-fabs(ms)-1&&mms==ms) { - + for (fi=0; fi>dipole_transition_x[k][ks][fi]; for (fi=0; fi>dipole_transition_y[k][ks][fi]; for (fi=0; fi>dipole_transition_z[k][ks][fi]; - + break; - + } else for (fi=0; fi>F>>F>>F; } } - - + + for(int fi=0;fifield_thresh[k]) {En[k] = 0; Gamman[k] = 0;} + if(Ez_stat>field_thresh[k]) {En[k] = 0; Gamman[k] = 0;} else { En[k] = energy[k][iEz]+cEz*(energy[k][iEz+1]-energy[k][iEz]); - + if(gamma_autoionization[k][iEz]<1e-21) Gamman[k] = 0; - else + else Gamman[k] = gamma_autoionization[k][iEz]*pow(gamma_autoionization[k][iEz+1]/gamma_autoionization[k][iEz],cEz); } } @@ -347,22 +347,22 @@ int Stark::getStates() {return st;} double Stark::getStarkEnergy( - double mass, int n1, int n2, int m, + double mass, int n1, int n2, int m, double Ex,double Ey,double Ez, double Bx,double By,double Bz, double px,double py,double pz){ -LorentzTransformationEM::transform(mass,px,py,pz,Ex,Ey,Ez,Bx,By,Bz); +//LorentzTransformationEM::transform(mass,px,py,pz,Ex,Ey,Ez,Bx,By,Bz); -SetE(sqrt(Ex*Ex+Ey*Ey+Ez*Ez)/5.14220642e011); +SetE(sqrt(Ex*Ex+Ey*Ey+Ez*Ez)/5.14220642e011); return En[convert3to1level(n1+n2+abs(m)+1,n1,m)]; } - - - - + + + + diff --git a/Stark/Stark.hh b/Stark/Stark.hh index 92629c3..be67025 100644 --- a/Stark/Stark.hh +++ b/Stark/Stark.hh @@ -2,8 +2,8 @@ #define STARK_HH_ -#include "tcomplex.hh" -#include "CppPyWrapper.hh" +#include +#include #include diff --git a/Stark/wrap_stark.cc b/Stark/wrap_stark.cc index cdf601e..4055e89 100644 --- a/Stark/wrap_stark.cc +++ b/Stark/wrap_stark.cc @@ -1,7 +1,6 @@ -#include "orbit_mpi.hh" -#include "pyORBIT_Object.hh" +#include +#include -#include "wrap_utils.hh" #include "wrap_stark.hh" #include @@ -12,7 +11,6 @@ #include "Stark.hh" //using namespace OrbitUtils; -using namespace wrap_orbit_utils; namespace wrap_stark{ @@ -45,28 +43,28 @@ extern "C" { if(!PyArg_ParseTuple( args,"si:",&addressEG,&states)){ error("Stark(address,states) - params. are needed"); - } + } else { - + std::string addr(addressEG); self->cpp_obj = new Stark(addr,states); ((Stark*) self->cpp_obj)->setPyWrapper((PyObject*) self); } - + return 0; - - + + } - - - - + + + + static PyObject* Stark_getStarkEnergy(PyObject *self, PyObject *args){ Stark* cpp_Stark = (Stark*)((pyORBIT_Object*) self)->cpp_obj; - - + + int nVars = PyTuple_Size(args); double val; double mass; @@ -82,7 +80,7 @@ extern "C" { double px; double py; double pz; - + //NO NEW OBJECT CREATED BY PyArg_ParseTuple! NO NEED OF Py_DECREF() if(!PyArg_ParseTuple( args,"diiiddddddddd:",&mass,&n1,&n2,&m,&E_x,&E_y,&E_z,&B_x,&B_y,&B_z,&px,&py,&pz)) @@ -91,12 +89,12 @@ extern "C" { val=cpp_Stark->getStarkEnergy(mass,n1,n2,m,E_x,E_y,E_z,B_x,B_y,B_z,px,py,pz); return Py_BuildValue("d",val); } - - - - - + + + + + //----------------------------------------------------- //destructor for python PyBaseFieldSource class (__del__ method). @@ -104,7 +102,7 @@ extern "C" { static void Stark_del(pyORBIT_Object* self){ //std::cerr<<"The CppBaseFieldSource __del__ has been called!"<cpp_obj); - self->ob_type->tp_free((PyObject*)self); + Py_TYPE(self)->tp_free((PyObject*)self); } // defenition of the methods of the python PyBaseFieldSource wrapper class @@ -130,8 +128,7 @@ extern "C" { //new python PyBaseFieldSource wrapper type definition static PyTypeObject pyORBIT_Stark_Type = { - PyObject_HEAD_INIT(NULL) - 0, /*ob_size*/ + PyVarObject_HEAD_INIT(NULL, 0) "Stark", /*tp_name*/ sizeof(pyORBIT_Object), /*tp_basicsize*/ 0, /*tp_itemsize*/ @@ -169,7 +166,7 @@ extern "C" { (initproc) Stark_init, /* tp_init */ 0, /* tp_alloc */ Stark_new, /* tp_new */ - }; + }; //-------------------------------------------------- //Initialization function of the pyPyBaseFieldSource class diff --git a/TwoLevelAtom/TwoLevelAtom.cc b/TwoLevelAtom/TwoLevelAtom.cc index aa9e311..30ee645 100644 --- a/TwoLevelAtom/TwoLevelAtom.cc +++ b/TwoLevelAtom/TwoLevelAtom.cc @@ -6,20 +6,20 @@ // CREATED // 06/27/2008 // -// DESCRIPTION -// +// DESCRIPTION +// // /////////////////////////////////////////////////////////////////////////// -#include "orbit_mpi.hh" +#include #include #include #include #include #include "TwoLevelAtom.hh" -#include "RungeKuttaTracker.hh" -#include "OrbitConst.hh" -#include "LorentzTransformationEM.hh" +#include +#include +//#include "LorentzTransformationEM.hh" #define Re1(i) AmplAttr->attArr(i)[1] //i-part index, n,m-attr index @@ -40,19 +40,19 @@ using namespace OrbitUtils; TwoLevelAtom::TwoLevelAtom(BaseLaserFieldSource* BaseLaserField, double delta_E, double dipole_tr) { setName("unnamed"); - + LaserField=BaseLaserField; d_Energy=delta_E; dip_transition=dipole_tr; - + if(LaserField->getPyWrapper() != NULL){ Py_INCREF(LaserField->getPyWrapper()); } } -TwoLevelAtom::~TwoLevelAtom() +TwoLevelAtom::~TwoLevelAtom() { - + if(LaserField->getPyWrapper() == NULL){ delete LaserField; } else { @@ -62,98 +62,98 @@ TwoLevelAtom::~TwoLevelAtom() } void TwoLevelAtom::CalcPopulations(int i, Bunch* bunch) { - + PopAttr->attArr(i)[1] = Re1(i)*Re1(i)+Im1(i)*Im1(i); PopAttr->attArr(i)[2] = Re2(i)*Re2(i)+Im2(i)*Im2(i); PopAttr->attArr(i)[0] = 1 - PopAttr->attArr(i)[1]-PopAttr->attArr(i)[2]; - + } void TwoLevelAtom::setupEffects(Bunch* bunch){ - - + + if(bunch->hasParticleAttributes("Amplitudes")==0) { std::map part_attr_dict; part_attr_dict["size"] = 5; bunch->addParticleAttributes("Amplitudes",part_attr_dict); - + for (int i=0; igetSize();i++){ bunch->getParticleAttributes("Amplitudes")->attValue(i,1) = 1; } } - - + + if(bunch->hasParticleAttributes("Populations")==0) { std::map part_attr_dict; part_attr_dict["size"] = 3; bunch->addParticleAttributes("Populations",part_attr_dict); } - - + + AmplAttr = bunch->getParticleAttributes("Amplitudes"); PopAttr = bunch->getParticleAttributes("Populations"); - + for (int i=0; igetSize();i++) CalcPopulations(i, bunch); } -void TwoLevelAtom::applyEffectsForEach(Bunch* bunch, int i, - double* y_in_vct, double* y_out_vct, - double t, double t_step, +void TwoLevelAtom::applyEffectsForEach(Bunch* bunch, int i, + double* y_in_vct, double* y_out_vct, + double t, double t_step, OrbitUtils::BaseFieldSource* fieldSource, RungeKuttaTracker* tracker) { - + x0 = y_in_vct[0];y0 = y_in_vct[1];z0 = y_in_vct[2]; px0 = y_in_vct[3];py0 = y_in_vct[4];pz0 = y_in_vct[5]; x = y_out_vct[0];y = y_out_vct[1];z = y_out_vct[2]; px = y_out_vct[3];py = y_out_vct[4];pz = y_out_vct[5]; - + if(LaserField->region(x,y,z)) { - // This function gives parameters Ez_stat Ex_las[1...3] Ey_las[1...3] Ez_las[1...3] - //in natural unts (Volt per meter) in the frame of particle + // This function gives parameters Ez_stat Ex_las[1...3] Ey_las[1...3] Ez_las[1...3] + //in natural unts (Volt per meter) in the frame of particle GetParticleFrameFields(i, t,t_step, bunch); - - // This function gives parameters Ez_stat Ex_las[1...3] Ey_las[1...3] Ez_las[1...3] t_part omega_part part_t_step - //in atomic units in frame of particle - GetParticleFrameParameters(i,t,t_step,bunch); - + + // This function gives parameters Ez_stat Ex_las[1...3] Ey_las[1...3] Ez_las[1...3] t_part omega_part part_t_step + //in atomic units in frame of particle + GetParticleFrameParameters(i,t,t_step,bunch); + // This function provides step solution for density matrix using rk4 method - AmplSolver4step(i,bunch); - + AmplSolver4step(i,bunch); + CalcPopulations(i, bunch); - + } - + } void TwoLevelAtom::GetParticleFrameFields(int i,double t,double t_step, Bunch* bunch) { double E_x,E_y,E_z,nx,ny,nz,Eabs; - + for (int j=0; j<3;j++) { - + exp_phasa[j] = LaserField->getLaserEMField(x0+j*(x-x0)/2,y0+j*(y-y0)/2,z0+j*(z-z0)/2, t+j*t_step/2,Ex_las[j],Ey_las[j],Ez_las[j],Bx_las[j],By_las[j],Bz_las[j]); - - LorentzTransformationEM::complex_electric_transform(bunch->getMass(), - px0,py0,pz0, - Ex_las[j],Ey_las[j],Ez_las[j], - Bx_las[j],By_las[j],Bz_las[j]); - + + //LorentzTransformationEM::complex_electric_transform(bunch->getMass(), + // px0,py0,pz0, + // Ex_las[j],Ey_las[j],Ez_las[j], + // Bx_las[j],By_las[j],Bz_las[j]); + if(j==0) { E_x=abs(Ex_las[j]); E_y=abs(Ey_las[j]); E_z=abs(Ez_las[j]); Eabs=sqrt(pow(E_x,2)+pow(E_y,2)+pow(E_z,2)); } - - Ez_las[j] = Eabs*exp_phasa[j]; + + Ez_las[j] = Eabs*exp_phasa[j]; } // cout<getMass(); - + //This line calculates relyativistic factor-Gamma double gamma=sqrt(m*m+px0*px0+py0*py0+pz0*pz0)/m; double coeff=1./(gamma*ta); part_t_step=t_step*coeff; //time step in frame of particle (in atomic units) - t_part=t*coeff; - + t_part=t*coeff; + //Convertion form SI units to atomic units - for(int j=0;j<3;j++) { - Ez_las[j]/=Ea; - + for(int j=0;j<3;j++) { + Ez_las[j]/=Ea; + } - //cout<getFrequencyOmega(m,x0(i),y0(i),z0(i),px0(i),py0(i),pz0(i),t)<<"\n"; + //cout<getFrequencyOmega(m,x0(i),y0(i),z0(i),px0(i),py0(i),pz0(i),t)<<"\n"; } void TwoLevelAtom::AmplSolver4step(int i, Bunch* bunch) { tcomplex z1,z2; double dt; - - for(int j=0; j<3;j++) mu_Elas[j]=dip_transition*conj(Ez_las[j])*J/2.; + + for(int j=0; j<3;j++) mu_Elas[j]=dip_transition*conj(Ez_las[j])*J/2.; z1=exp(J*t_part*fabs(d_Energy)); z2=exp(J*part_t_step*fabs(d_Energy)/2.); - exp_mu_El[1]=z1*mu_Elas[0]; + exp_mu_El[1]=z1*mu_Elas[0]; exp_mu_El[2]=z1*z2*mu_Elas[1]; exp_mu_El[3]=exp_mu_El[2]; exp_mu_El[4]=z1*z2*z2*mu_Elas[2]; - - for(int j=1; j<5; j++) { + + for(int j=1; j<5; j++) { if (j==4) dt=part_t_step; else dt=part_t_step/2.; - k_RungeKutt_1[j]=conj(exp_mu_El[j])*(Ampl_2(i)+k_RungeKutt_2[j-1]*dt); - k_RungeKutt_2[j]=-exp_mu_El[j]*(Ampl_1(i)+k_RungeKutt_1[j-1]*dt); + k_RungeKutt_1[j]=conj(exp_mu_El[j])*(Ampl_2(i)+k_RungeKutt_2[j-1]*dt); + k_RungeKutt_2[j]=-exp_mu_El[j]*(Ampl_1(i)+k_RungeKutt_1[j-1]*dt); } - - z1=(k_RungeKutt_1[1]+2.*k_RungeKutt_1[2]+2.*k_RungeKutt_1[3]+k_RungeKutt_1[4])/6.; + + z1=(k_RungeKutt_1[1]+2.*k_RungeKutt_1[2]+2.*k_RungeKutt_1[3]+k_RungeKutt_1[4])/6.; z2=(k_RungeKutt_2[1]+2.*k_RungeKutt_2[2]+2.*k_RungeKutt_2[3]+k_RungeKutt_2[4])/6.; - + Re1(i)+=z1.real()*part_t_step; Im1(i)+=z1.imag()*part_t_step; - + Re2(i)+=z2.real()*part_t_step; Im2(i)+=z2.imag()*part_t_step; } diff --git a/TwoLevelAtom/TwoLevelAtom.hh b/TwoLevelAtom/TwoLevelAtom.hh index 217fe40..3080efd 100644 --- a/TwoLevelAtom/TwoLevelAtom.hh +++ b/TwoLevelAtom/TwoLevelAtom.hh @@ -3,9 +3,9 @@ #include "Python.h" -#include "ExternalEffects.hh" +#include #include "BaseLaserFieldSource.hh" -#include "ParticleAttributes.hh" +#include using namespace TrackerRK4; diff --git a/TwoLevelAtom/wrap_two_level_atom.cc b/TwoLevelAtom/wrap_two_level_atom.cc index 6a1065c..2425c43 100644 --- a/TwoLevelAtom/wrap_two_level_atom.cc +++ b/TwoLevelAtom/wrap_two_level_atom.cc @@ -1,7 +1,6 @@ -#include "orbit_mpi.hh" -#include "pyORBIT_Object.hh" +#include +#include -#include "wrap_utils.hh" #include "wrap_two_level_atom.hh" #include @@ -47,7 +46,7 @@ extern "C" { if(!PyArg_ParseTuple( args,"Odd:",&pyBaseLaserField,&delta_E,&dipole_tr)){ error("TwoLevelAtom(LaserField,delta_E,dipole_transition) - params. are needed"); - } + } else { BaseLaserFieldSource* lfs = (BaseLaserFieldSource*) ((pyORBIT_Object*) pyBaseLaserField)->cpp_obj; @@ -57,7 +56,7 @@ extern "C" { return 0; } - // name([name]) - sets or returns the name of the External Effeects class + // name([name]) - sets or returns the name of the External Effeects class static PyObject* TwoLevelAtom_name(PyObject *self, PyObject *args){ TwoLevelAtom* cpp_TwoLevelAtom = (TwoLevelAtom*) ((pyORBIT_Object*) self)->cpp_obj; const char* name = NULL; @@ -69,15 +68,15 @@ extern "C" { cpp_TwoLevelAtom->setName(name_str); } return Py_BuildValue("s",cpp_TwoLevelAtom->getName().c_str()); - } - + } + //----------------------------------------------------- //destructor for python PyExternalEffects class (__del__ method). //----------------------------------------------------- static void TwoLevelAtom_del(pyORBIT_Object* self){ //std::cerr<<"The LasStripExternalEffects __del__ has been called!"<cpp_obj); - self->ob_type->tp_free((PyObject*)self); + Py_TYPE(self)->tp_free((PyObject*)self); } // defenition of the methods of the python PyExternalEffects wrapper class @@ -96,8 +95,7 @@ extern "C" { //new python PyExternalEffects wrapper type definition static PyTypeObject pyORBIT_TwoLevelAtom_Type = { - PyObject_HEAD_INIT(NULL) - 0, /*ob_size*/ + PyVarObject_HEAD_INIT(NULL, 0) "TwoLevelAtom", /*tp_name*/ sizeof(pyORBIT_Object), /*tp_basicsize*/ 0, /*tp_itemsize*/ @@ -135,10 +133,10 @@ extern "C" { (initproc) TwoLevelAtom_init, /* tp_init */ 0, /* tp_alloc */ TwoLevelAtom_new, /* tp_new */ - }; + }; + - //-------------------------------------------------- //Initialization function of the pyPyExternalEffects class //It will be called from Bunch wrapper initialization @@ -147,7 +145,7 @@ extern "C" { if (PyType_Ready(&pyORBIT_TwoLevelAtom_Type) < 0) return; Py_INCREF(&pyORBIT_TwoLevelAtom_Type); PyModule_AddObject(module, "TwoLevelAtom", (PyObject *)&pyORBIT_TwoLevelAtom_Type); - + } #ifdef __cplusplus diff --git a/TwoLevelStrongField/StarkStrongField.cc b/TwoLevelStrongField/StarkStrongField.cc index d0f1f66..2aded14 100644 --- a/TwoLevelStrongField/StarkStrongField.cc +++ b/TwoLevelStrongField/StarkStrongField.cc @@ -7,19 +7,19 @@ // 04/21/2003 // // DESCRIPTION -// The base class for Python implementation of a field source. -// It should be sub-classed on Python level and implements +// The base class for Python implementation of a field source. +// It should be sub-classed on Python level and implements // getElectricField(x,y,z,t) and getMagneticField (x,y,z,t) methods. // The results of these methods will be available from the c++ level. // This is an example of embedding Python in C++ Orbit level. // // /////////////////////////////////////////////////////////////////////////// -#include "orbit_mpi.hh" -#include "BufferStore.hh" +#include +#include #include "StarkStrongField.hh" -#include "LorentzTransformationEM.hh" +//#include "LorentzTransformationEM.hh" #include #include #include @@ -37,12 +37,12 @@ inline int convert3to1level(int n,int n1, int m){ StarkStrongField::StarkStrongField(std::string addressEG,int n1,int n2,int mm) -{ +{ int rank_MPI,size_MPI; ORBIT_MPI_Comm_size(MPI_COMM_WORLD, &size_MPI); ORBIT_MPI_Comm_rank(MPI_COMM_WORLD, &rank_MPI); - + std::ifstream file; std::ifstream file1; std::ifstream file2; @@ -59,22 +59,22 @@ StarkStrongField::StarkStrongField(std::string addressEG,int n1,int n2,int mm) - - - + + + if(rank_MPI == 0) { - - snprintf(nameEG,MAX_LENGTH_ADDRESS,"%sStarkEG/000.txt",addressEG.c_str()); - file.open(nameEG); file>>F>>F>>F>>delta_F; file.clear();file.close(); - snprintf(nameEG,MAX_LENGTH_ADDRESS,"%sStarkEG/%i%i%i.txt",addressEG.c_str(),n1,n2,abs(m)); - n_data=-1; file.open(nameEG); while(!file.eof()) {file>>field_thresh>>F>>G_ch; n_data++;} file.clear();file.close(); - + snprintf(nameEG,MAX_LENGTH_ADDRESS,"%sStarkEG/000.txt",addressEG.c_str()); + file.open(nameEG); file>>F>>F>>F>>delta_F; file.clear();file.close(); + + snprintf(nameEG,MAX_LENGTH_ADDRESS,"%sStarkEG/%i%i%i.txt",addressEG.c_str(),n1,n2,abs(m)); + n_data=-1; file.open(nameEG); while(!file.eof()) {file>>field_thresh>>F>>G_ch; n_data++;} file.clear();file.close(); - } + + } ORBIT_MPI_Bcast(&delta_F, 1, MPI_DOUBLE,0,MPI_COMM_WORLD); @@ -84,41 +84,41 @@ StarkStrongField::StarkStrongField(std::string addressEG,int n1,int n2,int mm) - - deltaE = new double[n_data]; - gamma_autoionization = new double[n_data]; + + deltaE = new double[n_data]; + gamma_autoionization = new double[n_data]; dipole_transition = new double[n_data]; - - - + + + if(rank_MPI == 0) { snprintf(nameEG,MAX_LENGTH_ADDRESS,"%sStarkEG/000.txt",addressEG.c_str()); - snprintf(nameEG1,MAX_LENGTH_ADDRESS,"%sStarkEG/%i%i%i.txt",addressEG.c_str(),n1,n2,abs(m)); + snprintf(nameEG1,MAX_LENGTH_ADDRESS,"%sStarkEG/%i%i%i.txt",addressEG.c_str(),n1,n2,abs(m)); snprintf(nameEG2,MAX_LENGTH_ADDRESS,"%sTDM_StrongField/TDM_%i%i%i.txt",addressEG.c_str(),n1,n2,abs(m)); file.open(nameEG); file1.open(nameEG1); file2.open(nameEG2); - for (fi=0; fi>F>>E0>>G_ch; + for (fi=0; fi>F>>E0>>G_ch; file1>>F>>E1>>G_ch; gamma_autoionization[fi] = atof(G_ch); file2>>F>>dipole_transition[fi]; - - deltaE[fi] = E1 - E0;} + + deltaE[fi] = E1 - E0;} file.clear();file.close(); file1.clear();file1.close(); file2.clear();file2.close(); } - - + + for(int fi=0;fifield_thresh) {deltaEn = 0; Gamman = 0; dipole_transition = 0;} + + if(E>field_thresh) {deltaEn = 0; Gamman = 0; dipole_transition = 0;} else { deltaEn = deltaE[iEz]+cEz*(deltaE[iEz+1]-deltaE[iEz]); dipole_transitionn = dipole_transition[iEz]+cEz*(dipole_transition[iEz+1]-dipole_transition[iEz]); - + if(gamma_autoionization[iEz]<1e-100) Gamman = 0; - else + else Gamman = gamma_autoionization[iEz]*pow(gamma_autoionization[iEz+1]/gamma_autoionization[iEz],cEz); } } - + double StarkStrongField::getStarkStrongFielddeltaE( - double mass, + double mass, double Ex,double Ey,double Ez, double Bx,double By,double Bz, double px,double py,double pz){ -LorentzTransformationEM::transform(mass,px,py,pz,Ex,Ey,Ez,Bx,By,Bz); +//LorentzTransformationEM::transform(mass,px,py,pz,Ex,Ey,Ez,Bx,By,Bz); -SetE(sqrt(Ex*Ex+Ey*Ey+Ez*Ez)/5.14220642e011); +SetE(sqrt(Ex*Ex+Ey*Ey+Ez*Ez)/5.14220642e011); return deltaEn; } - - - - + + + + diff --git a/TwoLevelStrongField/StarkStrongField.hh b/TwoLevelStrongField/StarkStrongField.hh index e95754b..0234099 100644 --- a/TwoLevelStrongField/StarkStrongField.hh +++ b/TwoLevelStrongField/StarkStrongField.hh @@ -2,8 +2,8 @@ #define STARKSTRONGFIELD_HH_ -#include "tcomplex.hh" -#include "CppPyWrapper.hh" +#include +#include #include diff --git a/TwoLevelStrongField/TwoLevelStrongField.cc b/TwoLevelStrongField/TwoLevelStrongField.cc index f19235e..a98da3b 100644 --- a/TwoLevelStrongField/TwoLevelStrongField.cc +++ b/TwoLevelStrongField/TwoLevelStrongField.cc @@ -6,20 +6,20 @@ // CREATED // 06/27/2008 // -// DESCRIPTION -// +// DESCRIPTION +// // /////////////////////////////////////////////////////////////////////////// -#include "orbit_mpi.hh" +#include #include #include #include #include #include "TwoLevelStrongField.hh" -#include "RungeKuttaTracker.hh" -#include "OrbitConst.hh" -#include "LorentzTransformationEM.hh" +#include +#include +//#include "LorentzTransformationEM.hh" #include "FieldRotation.hh" @@ -41,224 +41,224 @@ using namespace OrbitUtils; TwoLevelStrongField::TwoLevelStrongField(BaseLaserFieldSource* BaseLaserField, StarkStrongField* StarkStr) { setName("unnamed"); - + LaserField = BaseLaserField; Stark = StarkStr; - + if(LaserField->getPyWrapper() != NULL){ Py_INCREF(LaserField->getPyWrapper()); } - + if(Stark->getPyWrapper() != NULL){ Py_INCREF(Stark->getPyWrapper()); } - + } -TwoLevelStrongField::~TwoLevelStrongField() +TwoLevelStrongField::~TwoLevelStrongField() { - + if(LaserField->getPyWrapper() == NULL){ delete LaserField; } else { Py_XDECREF(LaserField->getPyWrapper()); } - + if(Stark->getPyWrapper() == NULL){ delete Stark; } else { Py_XDECREF(Stark->getPyWrapper()); } - - + + } void TwoLevelStrongField::CalcPopulations(int i, Bunch* bunch) { - + PopAttr->attArr(i)[1] = Re1(i)*Re1(i)+Im1(i)*Im1(i); PopAttr->attArr(i)[2] = Re2(i)*Re2(i)+Im2(i)*Im2(i); PopAttr->attArr(i)[0] = 1 - PopAttr->attArr(i)[1]-PopAttr->attArr(i)[2]; - + } void TwoLevelStrongField::setupEffects(Bunch* bunch){ - + int_E = new tcomplex[bunch->getSize()]; - + if(bunch->hasParticleAttributes("Amplitudes")==0) { std::map part_attr_dict; part_attr_dict["size"] = 5; bunch->addParticleAttributes("Amplitudes",part_attr_dict); - + for (int i=0; igetSize();i++){ bunch->getParticleAttributes("Amplitudes")->attValue(i,1) = 1; } } - - + + if(bunch->hasParticleAttributes("Populations")==0) { std::map part_attr_dict; part_attr_dict["size"] = 3; bunch->addParticleAttributes("Populations",part_attr_dict); } - - + + AmplAttr = bunch->getParticleAttributes("Amplitudes"); PopAttr = bunch->getParticleAttributes("Populations"); - - + + nx=new double[bunch->getSize()]; ny=new double[bunch->getSize()]; nz=new double[bunch->getSize()]; - + install_field_dir=new bool [bunch->getSize()]; - + for (int i=0; igetSize();i++) install_field_dir[i] = true; - - + + for (int i=0; igetSize();i++) CalcPopulations(i, bunch); } void TwoLevelStrongField::finalizeEffects(Bunch* bunch) { - + delete [] int_E; delete [] nx; delete [] ny; delete [] nz; - + delete [] install_field_dir; - + } -void TwoLevelStrongField::applyEffectsForEach(Bunch* bunch, int i, - double* y_in_vct, double* y_out_vct, - double t, double t_step, +void TwoLevelStrongField::applyEffectsForEach(Bunch* bunch, int i, + double* y_in_vct, double* y_out_vct, + double t, double t_step, OrbitUtils::BaseFieldSource* fieldSource, RungeKuttaTracker* tracker) { - + x0 = y_in_vct[0];y0 = y_in_vct[1];z0 = y_in_vct[2]; px0 = y_in_vct[3];py0 = y_in_vct[4];pz0 = y_in_vct[5]; x = y_out_vct[0];y = y_out_vct[1];z = y_out_vct[2]; px = y_out_vct[3];py = y_out_vct[4];pz = y_out_vct[5]; - - // This function gives parameters Ez_stat Ex_las[1...3] Ey_las[1...3] Ez_las[1...3] - //in natural unts (Volt per meter) in the frame of particle + + // This function gives parameters Ez_stat Ex_las[1...3] Ey_las[1...3] Ez_las[1...3] + //in natural unts (Volt per meter) in the frame of particle GetParticleFrameFields(i, t,t_step, bunch,fieldSource); - - // This function gives parameters Ez_stat Ex_las[1...3] Ey_las[1...3] Ez_las[1...3] t_part omega_part part_t_step - //in atomic units in frame of particle - GetParticleFrameParameters(i,t,t_step,bunch); - + + // This function gives parameters Ez_stat Ex_las[1...3] Ey_las[1...3] Ez_las[1...3] t_part omega_part part_t_step + //in atomic units in frame of particle + GetParticleFrameParameters(i,t,t_step,bunch); + // This function provides step solution for density matrix using rk4 method - AmplSolver4step(i,bunch); - + AmplSolver4step(i,bunch); + CalcPopulations(i, bunch); - + } void TwoLevelStrongField::GetParticleFrameFields(int i,double t,double t_step, Bunch* bunch, OrbitUtils::BaseFieldSource* fieldSource) { double E_x,E_y,E_z,Eabs; - - - fieldSource->getElectricMagneticField(x0,y0,z0,t,Ex_stat,Ey_stat,Ez_stat,Bx_stat,By_stat,Bz_stat); - LorentzTransformationEM::transform(bunch->getMass(),px0,py0,pz0,Ex_stat,Ey_stat,Ez_stat,Bx_stat,By_stat,Bz_stat); - + + + fieldSource->getElectricMagneticField(x0,y0,z0,t,Ex_stat,Ey_stat,Ez_stat,Bx_stat,By_stat,Bz_stat); + //LorentzTransformationEM::transform(bunch->getMass(),px0,py0,pz0,Ex_stat,Ey_stat,Ez_stat,Bx_stat,By_stat,Bz_stat); + if (install_field_dir[i]) { - nx[i] = Ex_stat; - ny[i] = Ey_stat; + nx[i] = Ex_stat; + ny[i] = Ey_stat; nz[i] = Ez_stat; install_field_dir[i]=false; - } - - + } + + Ez_stat=sqrt(Ex_stat*Ex_stat+Ey_stat*Ey_stat+Ez_stat*Ez_stat); - Ex_stat=0; + Ex_stat=0; Ey_stat=0; - - - + + + if(LaserField->region(x,y,z)) { - + for (int j=0; j<3;j++) { - + exp_phasa[j] = LaserField->getLaserEMField(x0+j*(x-x0)/2,y0+j*(y-y0)/2,z0+j*(z-z0)/2, t+j*t_step/2,Ex_las[j],Ey_las[j],Ez_las[j],Bx_las[j],By_las[j],Bz_las[j]); - - LorentzTransformationEM::complex_electric_transform(bunch->getMass(), - px0,py0,pz0, - Ex_las[j],Ey_las[j],Ez_las[j], - Bx_las[j],By_las[j],Bz_las[j]); - - + + //LorentzTransformationEM::complex_electric_transform(bunch->getMass(), + // px0,py0,pz0, + // Ex_las[j],Ey_las[j],Ez_las[j], + // Bx_las[j],By_las[j],Bz_las[j]); + + FieldRotation::RotateElectricFieldsV(nx[i],ny[i],nz[i],Ex_las[j],Ey_las[j],Ez_las[j]); - + if(j==0) { E_x=abs(Ex_las[j]); E_y=abs(Ey_las[j]); E_z=abs(Ez_las[j]); - - - + + + if(Stark->m == 0) Eabs = E_z; if(Stark->m == 1) Eabs=sqrt(pow(E_x,2)+pow(E_y,2)); } - - E_las[j] = Eabs*exp_phasa[j]; + + E_las[j] = Eabs*exp_phasa[j]; } } - else - for(int j=0;j<3;j++) { - E_las[j] = 0.; + else + for(int j=0;j<3;j++) { + E_las[j] = 0.; } - - + + // if(((int)(t/t_step))%1000 == 0) // cout<getMass(); - + //This line calculates relyativistic factor-Gamma double gamma=sqrt(m*m+px0*px0+py0*py0+pz0*pz0)/m; double coeff=1./(gamma*ta); part_t_step=t_step*coeff; //time step in frame of particle (in atomic units) - + //Convertion form SI units to atomic units - for(int j=0;j<3;j++) { - E_las[j]/=Ea; + for(int j=0;j<3;j++) { + E_las[j]/=Ea; } - - Ex_stat/=Ea; - Ey_stat/=Ea; + + Ex_stat/=Ea; + Ey_stat/=Ea; Ez_stat/=Ea; - + // if(((int)(t/t_step))%100 == 0) // cout<getFrequencyOmega(m,x0,y0,z0,px0,py0,pz0,t); - - + +// double omega = gamma*ta*LaserField->getFrequencyOmega(m,x0,y0,z0,px0,py0,pz0,t); + + } @@ -266,47 +266,47 @@ void TwoLevelStrongField::AmplSolver4step(int i, Bunch* bunch) { tcomplex z1,z2,z,JdEdt; double dt; - + Stark->SetE(Ez_stat); - + JdEdt = tcomplex(0.,Stark->deltaEn*part_t_step); - - - - - for(int j=0; j<3;j++) mu_Elas[j]=Stark->dipole_transitionn*conj(E_las[j])*J/2.; - - - - z1=exp(int_E[i]); + + + + + for(int j=0; j<3;j++) mu_Elas[j]=Stark->dipole_transitionn*conj(E_las[j])*J/2.; + + + + z1=exp(int_E[i]); z2=exp(JdEdt/2.); - exp_mu_El[1]=z1*mu_Elas[0]; + exp_mu_El[1]=z1*mu_Elas[0]; exp_mu_El[2]=z1*z2*mu_Elas[1]; exp_mu_El[3]=exp_mu_El[2]; exp_mu_El[4]=z1*z2*z2*mu_Elas[2]; - + int_E[i] += JdEdt; - - - for(int j=1; j<5; j++) { + + + for(int j=1; j<5; j++) { if (j==4) dt=part_t_step; else dt=part_t_step/2.; - k_RungeKutt_1[j]=conj(exp_mu_El[j])*(Ampl_2(i)+k_RungeKutt_2[j-1]*dt); - k_RungeKutt_2[j]=-exp_mu_El[j]*(Ampl_1(i)+k_RungeKutt_1[j-1]*dt); - + k_RungeKutt_1[j]=conj(exp_mu_El[j])*(Ampl_2(i)+k_RungeKutt_2[j-1]*dt); + k_RungeKutt_2[j]=-exp_mu_El[j]*(Ampl_1(i)+k_RungeKutt_1[j-1]*dt); + k_RungeKutt_2[j]-=Stark->Gamman*(Ampl_2(i)+k_RungeKutt_2[j-1]*dt)/2.; } - - - z1=(k_RungeKutt_1[1]+2.*k_RungeKutt_1[2]+2.*k_RungeKutt_1[3]+k_RungeKutt_1[4])/6.; + + + z1=(k_RungeKutt_1[1]+2.*k_RungeKutt_1[2]+2.*k_RungeKutt_1[3]+k_RungeKutt_1[4])/6.; z2=(k_RungeKutt_2[1]+2.*k_RungeKutt_2[2]+2.*k_RungeKutt_2[3]+k_RungeKutt_2[4])/6.; - + Re1(i)+=z1.real()*part_t_step; Im1(i)+=z1.imag()*part_t_step; - + Re2(i)+=z2.real()*part_t_step; Im2(i)+=z2.imag()*part_t_step; } diff --git a/TwoLevelStrongField/TwoLevelStrongField.hh b/TwoLevelStrongField/TwoLevelStrongField.hh index 4fdff20..c7fa2fa 100644 --- a/TwoLevelStrongField/TwoLevelStrongField.hh +++ b/TwoLevelStrongField/TwoLevelStrongField.hh @@ -3,9 +3,9 @@ #include "Python.h" -#include "ExternalEffects.hh" +#include #include "BaseLaserFieldSource.hh" -#include "ParticleAttributes.hh" +#include #include "StarkStrongField.hh" using namespace TrackerRK4; diff --git a/TwoLevelStrongField/wrap_stark_strong_field.cc b/TwoLevelStrongField/wrap_stark_strong_field.cc index 302d7ee..75aa836 100644 --- a/TwoLevelStrongField/wrap_stark_strong_field.cc +++ b/TwoLevelStrongField/wrap_stark_strong_field.cc @@ -1,7 +1,6 @@ -#include "orbit_mpi.hh" -#include "pyORBIT_Object.hh" +#include +#include -#include "wrap_utils.hh" #include "wrap_stark_strong_field.hh" #include @@ -12,7 +11,6 @@ #include "StarkStrongField.hh" //using namespace OrbitUtils; -using namespace wrap_orbit_utils; namespace wrap_stark_strong_field{ @@ -45,28 +43,28 @@ extern "C" { if(!PyArg_ParseTuple( args,"siii:",&addressEG,&n1,&n2,&m)){ error("StarkStrongField(address,n1,n2,m) - params. are needed"); - } + } else { - + std::string addr(addressEG); self->cpp_obj = new StarkStrongField(addr,n1,n2,m); ((StarkStrongField*) self->cpp_obj)->setPyWrapper((PyObject*) self); } - + return 0; - - + + } - - - - + + + + static PyObject* StarkStrongField_deltaE(PyObject *self, PyObject *args){ StarkStrongField* cpp_StarkStrongField = (StarkStrongField*)((pyORBIT_Object*) self)->cpp_obj; - - + + int nVars = PyTuple_Size(args); double val; double mass; @@ -79,7 +77,7 @@ extern "C" { double px; double py; double pz; - + //NO NEW OBJECT CREATED BY PyArg_ParseTuple! NO NEED OF Py_DECREF() if(!PyArg_ParseTuple( args,"dddddddddd:",&mass,&E_x,&E_y,&E_z,&B_x,&B_y,&B_z,&px,&py,&pz)) @@ -88,12 +86,12 @@ extern "C" { val=cpp_StarkStrongField->getStarkStrongFielddeltaE(mass,E_x,E_y,E_z,B_x,B_y,B_z,px,py,pz); return Py_BuildValue("d",val); } - - - - - + + + + + //----------------------------------------------------- //destructor for python PyBaseFieldSource class (__del__ method). @@ -101,7 +99,7 @@ extern "C" { static void StarkStrongField_del(pyORBIT_Object* self){ //std::cerr<<"The CppBaseFieldSource __del__ has been called!"<cpp_obj); - self->ob_type->tp_free((PyObject*)self); + Py_TYPE(self)->tp_free((PyObject*)self); } // defenition of the methods of the python PyBaseFieldSource wrapper class @@ -127,8 +125,7 @@ extern "C" { //new python PyBaseFieldSource wrapper type definition static PyTypeObject pyORBIT_StarkStrongField_Type = { - PyObject_HEAD_INIT(NULL) - 0, /*ob_size*/ + PyVarObject_HEAD_INIT(NULL, 0) "StarkStrongField", /*tp_name*/ sizeof(pyORBIT_Object), /*tp_basicsize*/ 0, /*tp_itemsize*/ @@ -166,7 +163,7 @@ extern "C" { (initproc) StarkStrongField_init, /* tp_init */ 0, /* tp_alloc */ StarkStrongField_new, /* tp_new */ - }; + }; //-------------------------------------------------- //Initialization function of the pyPyBaseFieldSource class diff --git a/TwoLevelStrongField/wrap_two_level_strong_field.cc b/TwoLevelStrongField/wrap_two_level_strong_field.cc index 8faf168..b4ffb7e 100644 --- a/TwoLevelStrongField/wrap_two_level_strong_field.cc +++ b/TwoLevelStrongField/wrap_two_level_strong_field.cc @@ -1,7 +1,6 @@ -#include "orbit_mpi.hh" -#include "pyORBIT_Object.hh" +#include +#include -#include "wrap_utils.hh" #include "wrap_two_level_strong_field.hh" #include @@ -48,7 +47,7 @@ extern "C" { if(!PyArg_ParseTuple( args,"OO:",&pyBaseLaserField,&pyStarkEffect)){ error("TwoLevelStrongField(LaserField,StarkEffect) - params. are needed"); - } + } else { BaseLaserFieldSource* lfs = (BaseLaserFieldSource*) ((pyORBIT_Object*) pyBaseLaserField)->cpp_obj; StarkStrongField* StarkEffect = (StarkStrongField*) ((pyORBIT_Object*) pyStarkEffect)->cpp_obj; @@ -59,7 +58,7 @@ extern "C" { return 0; } - // name([name]) - sets or returns the name of the External Effeects class + // name([name]) - sets or returns the name of the External Effeects class static PyObject* TwoLevelStrongField_name(PyObject *self, PyObject *args){ TwoLevelStrongField* cpp_TwoLevelStrongField = (TwoLevelStrongField*) ((pyORBIT_Object*) self)->cpp_obj; const char* name = NULL; @@ -71,15 +70,15 @@ extern "C" { cpp_TwoLevelStrongField->setName(name_str); } return Py_BuildValue("s",cpp_TwoLevelStrongField->getName().c_str()); - } - + } + //----------------------------------------------------- //destructor for python PyExternalEffects class (__del__ method). //----------------------------------------------------- static void TwoLevelStrongField_del(pyORBIT_Object* self){ //std::cerr<<"The LasStripExternalEffects __del__ has been called!"<cpp_obj); - self->ob_type->tp_free((PyObject*)self); + Py_TYPE(self)->tp_free((PyObject*)self); } // defenition of the methods of the python PyExternalEffects wrapper class @@ -98,8 +97,7 @@ extern "C" { //new python PyExternalEffects wrapper type definition static PyTypeObject pyORBIT_TwoLevelStrongField_Type = { - PyObject_HEAD_INIT(NULL) - 0, /*ob_size*/ + PyVarObject_HEAD_INIT(NULL, 0) "TwoLevelStrongField", /*tp_name*/ sizeof(pyORBIT_Object), /*tp_basicsize*/ 0, /*tp_itemsize*/ @@ -137,10 +135,10 @@ extern "C" { (initproc) TwoLevelStrongField_init, /* tp_init */ 0, /* tp_alloc */ TwoLevelStrongField_new, /* tp_new */ - }; + }; + - //-------------------------------------------------- //Initialization function of the pyPyExternalEffects class //It will be called from Bunch wrapper initialization @@ -149,7 +147,7 @@ extern "C" { if (PyType_Ready(&pyORBIT_TwoLevelStrongField_Type) < 0) return; Py_INCREF(&pyORBIT_TwoLevelStrongField_Type); PyModule_AddObject(module, "TwoLevelStrongField", (PyObject *)&pyORBIT_TwoLevelStrongField_Type); - + } #ifdef __cplusplus diff --git a/Walls/Walls.cc b/Walls/Walls.cc index ddbc756..1233194 100644 --- a/Walls/Walls.cc +++ b/Walls/Walls.cc @@ -7,26 +7,26 @@ // 06/27/2008 // // DESCRIPTION -// The base class for C++ implementation of a external effects -// during the transport of particles through the external field. +// The base class for C++ implementation of a external effects +// during the transport of particles through the external field. // It should be sub-classed on Python level and implement // setupEffects(Bunch* bunch) // finalizeEffects(Bunch* bunch) -// applyEffects(Bunch* bunch, int index, -// double* y_in_vct, double* y_out_vct, -// double t, double t_step, +// applyEffects(Bunch* bunch, int index, +// double* y_in_vct, double* y_out_vct, +// double t, double t_step, // OrbitUtils::BaseFieldSource* fieldSource) // methods. // The results of these methods will be available from the c++ level. // This is an example of embedding Python in C++ Orbit level. // -// -// -// -// +// +// +// +// // /////////////////////////////////////////////////////////////////////////// -#include "orbit_mpi.hh" +#include #include #include #include @@ -35,7 +35,7 @@ #include "Walls.hh" -#include "RungeKuttaTracker.hh" +#include @@ -71,7 +71,7 @@ Walls::Walls() -Walls::~Walls() +Walls::~Walls() { } @@ -92,9 +92,9 @@ Walls::~Walls() void Walls::setupEffects(Bunch* bunch){ - - + + if (bunch->hasParticleAttributes("pq_coords")==0) { std::map part_attr_dict; part_attr_dict["size"] = 6; @@ -109,23 +109,23 @@ void Walls::setupEffects(Bunch* bunch){ void Walls::prepareEffects(Bunch* bunch, double t){ - + for (int i=0; igetSize();i++) { x0(i)=bunch->coordArr()[i][0]; y0(i)=bunch->coordArr()[i][2]; z0(i)=bunch->coordArr()[i][4]; } - - + + } - - + + void Walls::finalizeEffects(Bunch* bunch) { - + if(bunch->hasParticleAttributes("pq_coords")==1) bunch->removeParticleAttributes("pq_coords"); - + } @@ -135,28 +135,28 @@ void Walls::finalizeEffects(Bunch* bunch) { /* bool Walls::crossPlane(int i,Bunch* bunch,double xs, double ys,double zs,double vx, double vy,double vz){ - + double vr = (x0(i)-x(i))*vx+(y0(i)-y(i))*vy+(z0(i)-z(i))*vz; - - + + if(vr==0) return false; - + else { double coeff = (vx*(xs-x0(i))+vy*(ys-y0(i))+vz*(zs-z0(i)))/vr; double x_cross = x0(i)+(x0(i)-x(i))*coeff; double y_cross = y0(i)+(y0(i)-y(i))*coeff; double z_cross = z0(i)+(z0(i)-z(i))*coeff; - + if ((x0(i)<=x_cross<=x(i))&&(y0(i)<=y_cross<=y(i))&&(z0(i)<=z_cross<=z(i))) return true; else return false; - - + + } - + } @@ -210,50 +210,50 @@ double R = 0.093115; cout<=x(i))&&(x(i)>=0.026719))) cout<<"foil"<<"\n"; - + if(((z0(i)<=z0)&&(z0<=z(i)) //intersection of flange z plane - + && - - (x1*x1+y1*y1<=R*R) //position of ring flange - + + (x1*x1+y1*y1<=R*R) //position of ring flange + && - + (x*x*gammax+2*ax*x*xp+bx*xp*xp < ex)&&(y*y*gammay+2*ay*y*yp+by*yp*yp < ey))) cout<<"ring"<<"\n"; */ if( - + ((z0(i)<=2.936809)&&(2.936809<=z(i))&&(0.118921>=x(i))&&(x(i)>=0.026719)) //position of thick foil - + || - + // ((z0(i)<=4.442887)&&(4.442887<=z(i))&&((x(i)+0.09679)*(x(i)+0.09679)+(y(i)-0)*(y(i)-0)<=0.0103476)) //position of ring flange - + ((z0(i)<=z0)&&(z0<=z(i)) //intersection of flange z plane - + && - - (x1*x1+y1*y1<=R*R) //position of ring flange - + + (x1*x1+y1*y1<=R*R) //position of ring flange + && - - (x*x*gammax+2*ax*x*xp+bx*xp*xp < ex)&&(y*y*gammay+2*ay*y*yp+by*yp*yp < ey)) //emittance of the ring beam - + + (x*x*gammax+2*ax*x*xp+bx*xp*xp < ex)&&(y*y*gammay+2*ay*y*yp+by*yp*yp < ey)) //emittance of the ring beam + // ((z0(i)<=2.936809)&&(2.936809<=z(i))&&(0.118921>=x(i))&&(x(i)>=0.026719)) - - + + // ((z0(i)<=0.25)&&(0.25<=z(i))) - + ) - - - + + + return true; else return false; @@ -266,7 +266,7 @@ else void Walls::applyEffects(Bunch* bunch, - double t, double t_step, + double t, double t_step, BaseFieldSource* fieldSource, RungeKuttaTracker* tracker) { @@ -276,9 +276,9 @@ void Walls::applyEffects(Bunch* bunch, if(crossSurface(i,bunch)) bunch->deleteParticleFast(i); - } + } + - } diff --git a/Walls/Walls.hh b/Walls/Walls.hh index 9ce79c2..fe0fe07 100644 --- a/Walls/Walls.hh +++ b/Walls/Walls.hh @@ -4,9 +4,9 @@ #include "Python.h" -#include "ExternalEffects.hh" +#include #include "BaseLaserFieldSource.hh" -#include "ParticleAttributes.hh" +#include diff --git a/Walls/wrap_walls.cc b/Walls/wrap_walls.cc index 4aca5d0..81aa35e 100644 --- a/Walls/wrap_walls.cc +++ b/Walls/wrap_walls.cc @@ -1,7 +1,6 @@ -#include "orbit_mpi.hh" -#include "pyORBIT_Object.hh" +#include +#include -#include "wrap_utils.hh" #include "wrap_walls.hh" #include @@ -39,24 +38,24 @@ extern "C" { //initializator for python PyExternalEffects class //this is implementation of the __init__ method static int Walls_init(pyORBIT_Object *self, PyObject *args, PyObject *kwds){ - - if(!PyArg_ParseTuple( args,"")){error("Walls() -No params. are needed");} - + + if(!PyArg_ParseTuple( args,"")){error("Walls() -No params. are needed");} + else { self->cpp_obj = new Walls(); ((Walls*) self->cpp_obj)->setPyWrapper((PyObject*) self); } - + return 0; } - - - - // name([name]) - sets or returns the name of the External Effeects class + + + + // name([name]) - sets or returns the name of the External Effeects class static PyObject* Walls_name(PyObject *self, PyObject *args){ Walls* cpp_Walls = (Walls*) ((pyORBIT_Object*) self)->cpp_obj; const char* name = NULL; @@ -68,18 +67,18 @@ extern "C" { cpp_Walls->setName(name_str); } return Py_BuildValue("s",cpp_Walls->getName().c_str()); - } - + } + + + - - //----------------------------------------------------- //destructor for python PyExternalEffects class (__del__ method). //----------------------------------------------------- static void Walls_del(pyORBIT_Object* self){ //std::cerr<<"The LasStripExternalEffects __del__ has been called!"<cpp_obj); - self->ob_type->tp_free((PyObject*)self); + Py_TYPE(self)->tp_free((PyObject*)self); } // defenition of the methods of the python PyExternalEffects wrapper class @@ -98,8 +97,7 @@ extern "C" { //new python PyExternalEffects wrapper type definition static PyTypeObject pyORBIT_Walls_Type = { - PyObject_HEAD_INIT(NULL) - 0, /*ob_size*/ + PyVarObject_HEAD_INIT(NULL, 0) "Walls", /*tp_name*/ sizeof(pyORBIT_Object), /*tp_basicsize*/ 0, /*tp_itemsize*/ @@ -137,10 +135,10 @@ extern "C" { (initproc) Walls_init, /* tp_init */ 0, /* tp_alloc */ Walls_new, /* tp_new */ - }; + }; + - //-------------------------------------------------- //Initialization function of the pyPyExternalEffects class //It will be called from Bunch wrapper initialization @@ -149,7 +147,7 @@ extern "C" { if (PyType_Ready(&pyORBIT_Walls_Type) < 0) return; Py_INCREF(&pyORBIT_Walls_Type); PyModule_AddObject(module, "Walls", (PyObject *)&pyORBIT_Walls_Type); - + } #ifdef __cplusplus diff --git a/examples/pyorbit3-1GeV-particle-stripping.py b/examples/pyorbit3-1GeV-particle-stripping.py new file mode 100644 index 0000000..153dbb1 --- /dev/null +++ b/examples/pyorbit3-1GeV-particle-stripping.py @@ -0,0 +1,120 @@ +#----------------------------------------------------- +# Track bunch with r and p through the external field +# The field is 1 T and has direction (0,1,0) +#----------------------------------------------------- + +import math + +from orbit.core.bunch import Bunch +from orbit.core.trackerrk4 import PyExternalEffects, RungeKuttaTracker, ExtEffectsContainer +from laserstripping import ConstEMfield +from laserstripping import HermiteGaussianLFmode, TwoLevelAtom + + +# ---------------------- Bunch/particle kinematic parameters ---------------------- # + +TK = 1.0 # energy of particle [GeV] +mass = 0.938256 + 0.000511 # mass [GeV] +charge = 0 # [e] +E = mass + TK +P = math.sqrt(E * E - mass * mass) +beta = P / E +# g = E / mass +c = 299792458 # speed of light +vz = c * beta # speed of particle + +b = Bunch() +b.charge(charge) +b.mass(mass) +b.addParticle(0, 0, 0, 0, 0, P) + + +# ---------------------- Bunch/particle quantum parameters ---------------------- # + +n_states = 3 # excitation between the 1st (ground) and the 3rd excited state +delta_E = 1. / 2. - 1. / (2. * n_states * n_states) # energy transition between 1st and 3rd excited state +dip_trans = math.sqrt( + 256 * n_states**7 * (n_states - 1)**(2 * n_states - 5) / + 3 / (n_states + 1)**(2 * n_states + 5) +) + +# ta = 2.418884326505e-17 # atomic unit of time + + +# ---------------------- Laser beam parameters ---------------------- # + +la0 = 2 * math.pi * 5.291772108e-11 / 7.297352570e-3 / delta_E # resonance wavelength at rest frame +la = 355.0e-9 # laser wavelength at the laboratory frame +te = TK - b.mass() * (la / la0 - 1) # temporal variable +kz = te / math.sqrt(P * P - te * te) # variable defining angle of interaction + +print("design angle of interaction =", math.atan2(1, -kz) * 360 / 2 / math.pi, + "direction of laser pulse =", (-1, 0, kz)) + +power = 1.0e6 # peak power of laser pulse in Watt +sigma_t = 25e-12 # rms time of Gaussian profile of laser pulse + +wx, wy, fx, fy = 0.001, 0.001, 1, 1 # default (non-critical) laser parameters +zl = -5 * sigma_t * c # initial shift of laser pulse before interaction + +LFS = HermiteGaussianLFmode( + math.sqrt(power), 0, 0, + abs(wx), abs(wy), fx, fy, + la, zl, sigma_t * c +) + +rx, ry, ax, ay = 0.2e-3, 0.2e-3, 0.5e-3, 0.5e-3 # rms size and divergence half angle at interaction point +LFS.setLocalParameters(abs(rx), abs(ry), ax, ay) + +# Orientation of beam (k-vector) {-1, 0, kz} and E-field polarization {kz, 0, 1} +LFS.setLaserFieldOrientation( + 0., 0., 0., + -1., 0., kz, + kz, 0., 1., + kz, 0., 1. +) +# LFS.setLaserFieldOrientation(0., 0., 0., -1., 0., kz, kz, 0., 1., 0., 1., 0.) # parallel polarization + +eff = TwoLevelAtom(LFS, delta_E, dip_trans) # excitation effect (ground state by default) + + +class ExternalEffects(PyExternalEffects): + def __init__(self): + PyExternalEffects.__init__(self) + + def setupEffects(self, bunch): + pass + + def prepareEffects(self, bunch, t): + pass + + def finalizeEffects(self, bunch): + pass + + def applyEffects(self, bunch, t, t_step, field_source, tracker): + if int(t / t_step) % (n_step // 1000) == 0: + print(b.z(0), "\t", b.partAttrValue("Populations", 0, 2)) + + def applyEffectsForEach(self, bunch, index, inVct, outVct, t, t_step, field_source, tracker): + pass + + +print_excitation_track = ExternalEffects() +cont_eff = ExtEffectsContainer() +# cont_eff.AddEffect(print_excitation_track) # prints excitation as a function of z +cont_eff.AddEffect(eff) + +b.z(0, b.z(0) - 5 * sigma_t * vz) # initial shift (z-coordinate) before interaction + +# Tracking time is supposed to be 2*(5*sigma_t) +n_step = 1000000 # number of tracking steps +time_step = 2 * (5 * sigma_t) / n_step # time step of tracking + +tracker = RungeKuttaTracker(0) +fS = ConstEMfield(0, 0, 0, 0, 0, 0) # permanent field E=0, B=0 (formally required by tracker) +tracker.track(b, 0, time_step * n_step, time_step, fS, cont_eff) # track from 0 for time_step*n_step + +print(b.partAttrValue("Populations", 0, 1)) # ground state population after excitation +print(b.partAttrValue("Populations", 0, 2)) # 3rd excited state after excitation +print(b.partAttrValue("Populations", 0, 1) + + b.partAttrValue("Populations", 0, 2)) # sum must equal 1 for TwoLevelAtom model diff --git a/meson.build b/meson.build new file mode 100644 index 0000000..ae6000d --- /dev/null +++ b/meson.build @@ -0,0 +1,102 @@ +project('laser-stripping', 'cpp') + +python = import('python').find_installation('python3') + +pyorbit3_dep = dependency('pyorbit3', required: true) +pyorbit3_incbase = pyorbit3_dep.get_variable(pkgconfig: 'includedir') / 'pyorbit3' + +# pyorbit3 headers reference each other with bare filenames like "Bunch.hh"; +# add all its subdirectories so those transitive includes resolve. +pyorbit3_subdirs_result = run_command( + 'find', pyorbit3_incbase, '-type', 'd', + check: false +).stdout().strip().split('\n') + +pyorbit3_subdir_incs = [] +foreach d : pyorbit3_subdirs_result + pyorbit3_subdir_incs += include_directories(d) +endforeach + +ls_sources = files([ + 'ConeLFMode/ConeLFmode.cc', + 'ConeLFMode/wrap_cone_lf_mode.cc', + 'ContinuumSS/ContinuumSS.cc', + 'ContinuumSS/TDMcontinuum.cc', + 'ContinuumSS/wrap_continuum_ss.cc', + 'ContinuumSS/wrap_tdm_continuum.cc', + 'DensityMatrix/DensityMatrix.cc', + 'DensityMatrix/wrap_density_matrix.cc', + 'DiagonalDensityMatrix/DM_noLaserField.cc', + 'DiagonalDensityMatrix/wrap_DM_noLaserField.cc', + 'FringeField/FringeField.cc', + 'FringeField/wrap_fringe_field.cc', + 'FroissartStoraLF/FroissartStoraLF.cc', + 'FroissartStoraLF/wrap_froissart_stora_lf.cc', + 'HermiteGaussianLFmode/HermiteGaussianLFmode.cc', + 'HermiteGaussianLFmode/wrap_hermite_gaussian_lf_mode.cc', + 'HminusStripping/HminusStripping.cc', + 'HminusStripping/wrap_hminus_stripping.cc', + 'HydrogenStarkParam/HydrogenStarkParam.cc', + 'HydrogenStarkParam/wrap_hydrogen_stark_param.cc', + 'LS_utils/BaseLaserFieldSource.cc', + 'LS_utils/FieldOrientation.cc', + 'LS_utils/FieldRotation.cc', + 'LS_utils/wrap_laserstripping.cc', + 'LaserFieldContainer/LaserFieldContainer.cc', + 'LaserFieldContainer/wrap_las_field_container.cc', + 'Quad_EM_field/ConstEMfield.cc', + 'Quad_EM_field/QuadEMfield.cc', + 'Quad_EM_field/wrap_const_em_field.cc', + 'Quad_EM_field/wrap_quad_em_field.cc', + 'Record/PrintExtEffects.cc', + 'Record/RecordEvolution.cc', + 'Record/wrap_print_ext_effects.cc', + 'Record/wrap_record_evolution.cc', + 'RegularGridFS/RegularGridFS.cc', + 'RegularGridFS/wrap_regular_grid_fs.cc', + 'SNSstrippingMagnet/SNSstrippingMagnet.cc', + 'SNSstrippingMagnet/wrap_snsstrippingmagnet.cc', + 'SchrodingerEquation/SchrodingerEquation.cc', + 'SchrodingerEquation/wrap_schrodinger_equation.cc', + 'Stark/Stark.cc', + 'Stark/wrap_stark.cc', + 'TwoLevelAtom/TwoLevelAtom.cc', + 'TwoLevelAtom/wrap_two_level_atom.cc', + 'TwoLevelStrongField/StarkStrongField.cc', + 'TwoLevelStrongField/TwoLevelStrongField.cc', + 'TwoLevelStrongField/wrap_stark_strong_field.cc', + 'TwoLevelStrongField/wrap_two_level_strong_field.cc', + 'Walls/Walls.cc', + 'Walls/wrap_walls.cc', +]) + +ls_inc = include_directories([ + '.', + 'ConeLFMode', + 'ContinuumSS', + 'DensityMatrix', + 'DiagonalDensityMatrix', + 'FringeField', + 'FroissartStoraLF', + 'HermiteGaussianLFmode', + 'HminusStripping', + 'HydrogenStarkParam', + 'LS_utils', + 'LaserFieldContainer', + 'Quad_EM_field', + 'Record', + 'RegularGridFS', + 'SNSstrippingMagnet', + 'SchrodingerEquation', + 'Stark', + 'TwoLevelAtom', + 'TwoLevelStrongField', + 'Walls', +]) + +python.extension_module('laserstripping', + sources: ls_sources, + include_directories: [ls_inc, pyorbit3_subdir_incs], + dependencies: [pyorbit3_dep], + install: true, +)