diff --git a/py/orbit/py_linac/lattice/LinacAccNodes.py b/py/orbit/py_linac/lattice/LinacAccNodes.py index d1119d8f..edeadb8f 100755 --- a/py/orbit/py_linac/lattice/LinacAccNodes.py +++ b/py/orbit/py_linac/lattice/LinacAccNodes.py @@ -297,13 +297,15 @@ def __init__(self, name = "quad"): self.setType("linacQuad") def fringeIN(node,paramsDict): - # B*rho = 3.335640952*momentum [T*m] if momentum in GeV/c + # B*rho = 3.335640952*momentum/charge [T*m] if momentum in GeV/c usageIN = node.getUsage() if(not usageIN): return - bunch = paramsDict["bunch"] + bunch = paramsDict["bunch"] + charge = bunch.charge() momentum = bunch.getSyncParticle().momentum() - kq = node.getParam("dB/dr")/(3.335640952*momentum) + #---- The charge sign will be accounted for inside tracking module functions. + kq = node.getParam("dB/dr")/bunch.B_Rho() poleArr = node.getParam("poles") klArr = node.getParam("kls") skewArr = node.getParam("skews") @@ -318,12 +320,15 @@ def fringeIN(node,paramsDict): TPB.multpfringeIN(bunch,pole,k,skew) def fringeOUT(node,paramsDict): + # B*rho = 3.335640952*momentum/charge [T*m] if momentum in GeV/c usageOUT = node.getUsage() if(not usageOUT): return bunch = paramsDict["bunch"] + charge = bunch.charge() momentum = bunch.getSyncParticle().momentum() - kq = node.getParam("dB/dr")/(3.335640952*momentum) + #---- The charge sign will be accounted for inside tracking module functions + kq = node.getParam("dB/dr")/bunch.B_Rho() poleArr = node.getParam("poles") klArr = node.getParam("kls") skewArr = node.getParam("skews") @@ -400,9 +405,14 @@ def track(self, paramsDict): The Quad Combined Function TEAPOT class implementation of the AccNode class track(probe) method. """ - bunch = paramsDict["bunch"] + bunch = paramsDict["bunch"] + charge = bunch.charge() momentum = bunch.getSyncParticle().momentum() - kq = self.getParam("dB/dr")/(3.335640952*momentum) + #---- The sign of dB/dr will be delivered to tracking module + #---- functions as kq. + #---- The charge sign will be accounted for inside tracking module + #---- functions. + kq = self.getParam("dB/dr")/bunch.B_Rho() nParts = self.getnParts() index = self.getActivePartIndex() length = self.getLength(index) @@ -646,8 +656,7 @@ def track(self, paramsDict): TPB.bend2(bunch, length) TPB.bend1(bunch, length, theta/2.0) return - - + class DCorrectorH(LinacMagnetNode): """ The Horizontal Dipole Corrector. @@ -684,11 +693,11 @@ def track(self, paramsDict): length = self.getParam("effLength")/nParts field = self.getParam("B") bunch = paramsDict["bunch"] + charge = bunch.charge() syncPart = bunch.getSyncParticle() momentum = syncPart.momentum() # dp/p = Q*c*B*L/p p in GeV/c c = 2.99792*10^8/10^9 - # Q is used inside kick-method - kick = field*length*0.299792/momentum + kick = -field*charge*length*0.299792/momentum self.tracking_module.kick(bunch,kick,0.,0.) class DCorrectorV(LinacMagnetNode): @@ -727,11 +736,11 @@ def track(self, paramsDict): length = self.getParam("effLength")/nParts field = self.getParam("B") bunch = paramsDict["bunch"] + charge = bunch.charge() syncPart = bunch.getSyncParticle() momentum = syncPart.momentum() # dp/p = Q*c*B*L/p p in GeV/c, c = 2.99792*10^8/10^9 - # Q is used inside kick-method - kick = field*length*0.299792/momentum + kick = field*charge*length*0.299792/momentum self.tracking_module.kick(bunch,0,kick,0.) class ThickKick(LinacMagnetNode): @@ -785,7 +794,8 @@ def track(self, paramsDict): """ The Thick Kick class implementation of the AccNode class track(probe) method. """ - bunch = paramsDict["bunch"] + bunch = paramsDict["bunch"] + charge = bunch.charge() momentum = bunch.getSyncParticle().momentum() Bx = self.getParam("Bx") By = self.getParam("By") @@ -795,9 +805,8 @@ def track(self, paramsDict): #print "debug name =",self.getName()," Bx=",Bx," By=",By," L=",self.getLength(index)," index=",index #========================================== # dp/p = Q*c*B*L/p p in GeV/c, c = 2.99792*10^8/10^9 - # Q is used inside kick-method - kickY = Bx*length*0.299792/momentum - kickX = By*length*0.299792/momentum + kickY = +Bx*charge*length*0.299792/momentum + kickX = -By*charge*length*0.299792/momentum self.tracking_module.drift(bunch, length/2.0) self.tracking_module.kick(bunch,kickX,kickY,0.) self.tracking_module.drift(bunch, length/2.0) diff --git a/py/orbit/py_linac/lattice/LinacFieldOverlappingNodes.py b/py/orbit/py_linac/lattice/LinacFieldOverlappingNodes.py index 49159b7f..63435c9f 100755 --- a/py/orbit/py_linac/lattice/LinacFieldOverlappingNodes.py +++ b/py/orbit/py_linac/lattice/LinacFieldOverlappingNodes.py @@ -283,7 +283,8 @@ def track(self, paramsDict): nParts = self.getnParts() index = self.getActivePartIndex() part_length = self.getLength(index) - bunch = paramsDict["bunch"] + bunch = paramsDict["bunch"] + charge = bunch.charge() syncPart = bunch.getSyncParticle() eKin_in = syncPart.kinEnergy() momentum = syncPart.momentum() @@ -308,10 +309,10 @@ def track(self, paramsDict): Ep = self.getEzFiledInternal(zp,rfCavity,E0L,rf_ampl) #------- track through a quad G = self.getTotalField((zm+z0)/2) - GP = 0. - if(self.useLongField == True): GP = self.getTotalFieldDerivative((zm+z0)/2) + dB_dz = 0. + if(self.useLongField == True): dB_dz = self.getTotalFieldDerivative((zm+z0)/2) if(abs(G) != 0.): - kq = G/(3.335640952*momentum) + kq = G/bunch.B_Rho() #------- track through a quad step = part_length/2 self.tracking_module.quad1(bunch,step/4.0, kq) @@ -319,9 +320,8 @@ def track(self, paramsDict): self.tracking_module.quad1(bunch,step/2.0, kq) self.tracking_module.quad2(bunch,step/2.0) self.tracking_module.quad1(bunch,step/4.0, kq) - if(abs(GP) != 0.): - kqP = GP/(3.335640952*momentum) - self.tracking_module.quad3(bunch,step, kqP) + if(abs(dB_dz) != 0.): + self.tracking_module.quad3(bunch,step,dB_dz) else: self.tracking_module.drift(bunch,part_length/2) self.part_pos += part_length/2 @@ -340,19 +340,18 @@ def track(self, paramsDict): self.cppGapModel.trackBunch(bunch,part_length/2,Em,E0,Ep,frequency,phase+delta_phase+modePhase) #------- track through a quad G = self.getTotalField((z0+zp)/2) - GP = 0. - if(self.useLongField == True): GP = self.getTotalFieldDerivative((z0+zp)/2) + dB_dz = 0. + if(self.useLongField == True): dB_dz = self.getTotalFieldDerivative((z0+zp)/2) if(abs(G) != 0.): - kq = G/(3.335640952*momentum) + kq = G/bunch.B_Rho() step = part_length/2 self.tracking_module.quad1(bunch,step/4.0, kq) self.tracking_module.quad2(bunch,step/2.0) self.tracking_module.quad1(bunch,step/2.0, kq) self.tracking_module.quad2(bunch,step/2.0) self.tracking_module.quad1(bunch,step/4.0, kq) - if(abs(GP) != 0.): - kqP = GP/(3.335640952*momentum) - self.tracking_module.quad3(bunch,step, kqP) + if(abs(dB_dz) != 0.): + self.tracking_module.quad3(bunch,step,dB_dz) else: self.tracking_module.drift(bunch,part_length/2) #---- advance the particle position @@ -615,15 +614,16 @@ def track(self, paramsDict): length = self.getLength(index) if(index == 0): self.z_value = - self.getLength()/2 bunch = paramsDict["bunch"] + charge = bunch.charge() momentum = bunch.getSyncParticle().momentum() n_steps = int(length/self.z_step)+1 z_step = length/n_steps for z_ind in range(n_steps): z = self.z_value + z_step*(z_ind+0.5) G = self.getTotalField(z) - GP = 0. - if(self.useLongField == True): GP = self.getTotalFieldDerivative(z) - kq = G/(3.335640952*momentum) + dB_dz = 0. + if(self.useLongField == True): dB_dz = self.getTotalFieldDerivative(z) + kq = G/bunch.B_Rho() if(abs(kq) == 0.): self.tracking_module.drift(bunch,z_step) continue @@ -633,9 +633,8 @@ def track(self, paramsDict): self.tracking_module.quad1(bunch,z_step/2.0, kq) self.tracking_module.quad2(bunch,z_step/2.0) self.tracking_module.quad1(bunch,z_step/4.0, kq) - if(abs(GP) != 0.): - kqP = GP/(3.335640952*momentum) - self.tracking_module.quad3(bunch,z_step, kqP) + if(abs(dB_dz) != 0.): + self.tracking_module.quad3(bunch,z_step, dB_dz) self.z_value += length def getTotalField(self,z_from_center): diff --git a/src/linac/tracking/linac_tracking.cc b/src/linac/tracking/linac_tracking.cc index 5de037d7..3ed38749 100644 --- a/src/linac/tracking/linac_tracking.cc +++ b/src/linac/tracking/linac_tracking.cc @@ -115,7 +115,7 @@ namespace linac_tracking // bunch = reference to the macro-particle bunch // length = length of transport // kq = quadrupole field strength [m^(-2)] - // kq = G/(B*rho) B*rho = 3.33564*P0 [T*m] where P0 in GeV/c + // kq = G/(B*rho) B*rho = 3.33564*P0/charge [T*m] where P0 in GeV/c // G = dB/dr [T/m] quad gradient // // RETURNS @@ -125,11 +125,7 @@ namespace linac_tracking void linac_quad1(Bunch* bunch, double length, double kq, int useCharge) { - double charge = +1.0; - if(useCharge == 1) charge = bunch->getCharge(); - - double kqc = kq * charge; - if(kqc == 0.) + if(kq == 0. || bunch->getCharge() == 0.) { linac_drift(bunch,length); return; @@ -154,8 +150,8 @@ namespace linac_tracking syncPart->setTime(syncPart->getTime() + delta_t); } - // ==== B*rho = 3.335640952*momentum [T*m] if momentum in GeV/c === - double dB_dr = kqc*3.335640952*syncPart->getMomentum(); + // ==== B*rho = 3.335640952*momentum/charge [T*m] if momentum in GeV/c === + double dB_dr = kq*bunch->getB_Rho(); double mass = syncPart->getMass(); double Ekin_s = syncPart->getEnergy(); @@ -175,6 +171,8 @@ namespace linac_tracking double** arr = bunch->coordArr(); int nParts = bunch->getSize(); + double kqc = 0.; + for(int i = 0; i < nParts; i++) { @@ -183,7 +181,7 @@ namespace linac_tracking Etotal = Ekin + mass; p2 = Ekin*(Ekin+2.0*mass); p = sqrt(p2); - kqc = dB_dr/(3.335640952*p); + kqc = dB_dr/(3.335640952*p/bunch->getCharge()); beta_z = p/Etotal; sqrt_kq = sqrt(fabs(kqc)); @@ -263,34 +261,34 @@ namespace linac_tracking // // DESCRIPTION // Quadrupole element 3: longitudinal field component of the quad field - // Bz(z) = x*y*dG(z)/dz + // Bz(x,y,z) = x*y*dG(z)/dz // This function performs the transverse kicks correction, so the length // is just a parameter. // // PARAMETERS // bunch = reference to the macro-particle bunch // length = length of transport - // kq = quadrupole field strength derivative [m^(-2)] - // kq = (dG/dz)/(B*rho) B*rho = 3.33564*P0 [T*m] where P0 in GeV/c // G = dB/dr [T/m] quad gradient + // 0.299792458 = (1/v_light) * 10^9 // // RETURNS // Nothing // /////////////////////////////////////////////////////////////////////////// - void linac_quad3(Bunch* bunch, double length, double kq, int useCharge) + void linac_quad3(Bunch* bunch, double length, double dB_dz) { - double charge = +1.0; - if(useCharge == 1) charge = bunch->getCharge(); - - double kqc = kq * charge; - if(kqc == 0.) + double charge = bunch->getCharge(); + + if(dB_dz == 0. || charge == 0.) { return; } - - double kick_coeff = kqc*length; + + SyncPart* syncPart = bunch->getSyncPart(); + //momentum in GeV/c + double momentum = syncPart->getMomentum(); + double kick_coeff = 0.299792458*charge*dB_dz*length/momentum; //coordinate array [part. index][x,xp,y,yp,z,dE] double** arr = bunch->coordArr(); @@ -325,12 +323,6 @@ namespace linac_tracking void kick(Bunch* bunch, double kx, double ky, double kE, int useCharge) { - double charge = +1.0; - if(useCharge == 1) charge = bunch->getCharge(); - - double kxc = kx * charge; - double kyc = ky * charge; - SyncPart* syncPart = bunch->getSyncPart(); double mass = syncPart->getMass(); @@ -345,7 +337,7 @@ namespace linac_tracking //coordinate array [part. index][x,xp,y,yp,z,dE] double** arr = bunch->coordArr(); - if(kxc != 0.) + if(kx != 0.) { for(int i = 0; i < bunch->getSize(); i++) { @@ -354,10 +346,10 @@ namespace linac_tracking p2 = Ekin*(Ekin+2.0*mass); p_z = sqrt(p2 - (arr[i][1]*arr[i][1] + arr[i][3]*arr[i][3])*p2_s); coeff = p_s/p_z; - arr[i][1] += kxc*coeff; + arr[i][1] += kx*coeff; } } - if(kyc != 0.) + if(ky != 0.) { for(int i = 0; i < bunch->getSize(); i++) { @@ -366,7 +358,7 @@ namespace linac_tracking p2 = Ekin*(Ekin+2.0*mass); p_z = sqrt(p2 - (arr[i][1]*arr[i][1] + arr[i][3]*arr[i][3])*p2_s); coeff = p_s/p_z; - arr[i][3] += kyc*coeff; + arr[i][3] += ky*coeff; } } if(kE != 0.) diff --git a/src/linac/tracking/linac_tracking.hh b/src/linac/tracking/linac_tracking.hh index aba2d24d..415225e9 100644 --- a/src/linac/tracking/linac_tracking.hh +++ b/src/linac/tracking/linac_tracking.hh @@ -46,7 +46,7 @@ namespace linac_tracking This function performs the transverse kicks correction, so the length is just a parameter. */ - void linac_quad3(Bunch* bunch, double length, double kq, int useCharge); + void linac_quad3(Bunch* bunch, double length, double dB_dz); /** diff --git a/src/linac/tracking/wrap_linac_tracking.cc b/src/linac/tracking/wrap_linac_tracking.cc index 2a17c93e..b6bb255d 100644 --- a/src/linac/tracking/wrap_linac_tracking.cc +++ b/src/linac/tracking/wrap_linac_tracking.cc @@ -75,15 +75,15 @@ extern "C" { PyObject* pyBunch; double length; - double kq; + double dB_dz = 0.; int useCharge = 1; if(!PyArg_ParseTuple( args, "Odd:linac_quad3", - &pyBunch, &length,&kq)) + &pyBunch, &length,&dB_dz)) { - error("linac tracking - linac_quad3(pyBunch,length,kq) - cannot parse arguments!"); + error("linac tracking - linac_quad3(pyBunch,length,dB_dz) - cannot parse arguments!"); } Bunch* cpp_bunch = (Bunch*) ((pyORBIT_Object *) pyBunch)->cpp_obj; - linac_tracking::linac_quad3(cpp_bunch, length, kq, useCharge); + linac_tracking::linac_quad3(cpp_bunch, length, dB_dz); Py_INCREF(Py_None); return Py_None; } diff --git a/src/orbit/Bunch.cc b/src/orbit/Bunch.cc index 6f635f12..f97ccef2 100644 --- a/src/orbit/Bunch.cc +++ b/src/orbit/Bunch.cc @@ -771,7 +771,20 @@ void Bunch::compress() double Bunch::getMass(){ return mass;} double Bunch::getCharge(){ return charge;} double Bunch::getClassicalRadius(){ return classicalRadius;} -double Bunch::getMacroSize(){ return macroSizeForAll;} +double Bunch::getMacroSize(){ return macroSizeForAll;} + +/** + Return the B*Rho parameter of the particle + B*Rho = momentum/charge - [T*m] + or B*Rho = 3.335640952*momentum[GeV/c]/charge[electron charges] +*/ +double Bunch::getB_Rho(){ + double B_Rho = 1.0e+36; + if(charge != 0. && syncPart != NULL){ + B_Rho = 3.335640952*syncPart->getMomentum()/charge; + } + return B_Rho; +} double Bunch::setMass(double val){ mass = val; diff --git a/src/orbit/Bunch.hh b/src/orbit/Bunch.hh index 2db1bba3..2d49e802 100644 --- a/src/orbit/Bunch.hh +++ b/src/orbit/Bunch.hh @@ -111,6 +111,13 @@ public: double getClassicalRadius(); // m double getCharge(); // sign and value in abs(e-charge) only double getMacroSize(); + + /** + Return the B*Rho parameter of the particle + B*Rho = momentum/charge - [T*m] + or B*Rho = 3.335640952*momentum[GeV/c]/charge[electron charges] + */ + double getB_Rho(); double setMass(double mass); // GeV double setCharge(double chrg); // sign and value in abs(e-charge) only diff --git a/src/orbit/wrap_bunch.cc b/src/orbit/wrap_bunch.cc index 8392e59b..4e5d4f2d 100644 --- a/src/orbit/wrap_bunch.cc +++ b/src/orbit/wrap_bunch.cc @@ -504,6 +504,13 @@ namespace wrap_orbit_bunch{ double val = cpp_bunch->getClassicalRadius(); return Py_BuildValue("d",val); } + + //Returns B_Rho of the particle in [Tesla*meter]. Parameter is used in TEAPOT + static PyObject* Bunch_B_Rho(PyObject *self, PyObject *args){ + Bunch* cpp_bunch = (Bunch*) ((pyORBIT_Object *) self)->cpp_obj; + double val = cpp_bunch->getB_Rho(); + return Py_BuildValue("d",val); + } //Sets or returns charge of the macro-particle in e-charge // the action is depended on the number of arguments @@ -1189,6 +1196,7 @@ namespace wrap_orbit_bunch{ { "ringwrap", Bunch_ringwrap ,METH_VARARGS,"Perform the ring wrap. Usage: ringwrap(ring_length)"}, { "mass", Bunch_mass ,METH_VARARGS,"Set mass(value) or get mass() the mass of particle in MeV"}, { "classicalRadius", Bunch_classicalRadius ,METH_VARARGS,"Returns a classical radius of particle in [m]"}, + { "B_Rho", Bunch_B_Rho ,METH_VARARGS,"Returns B*Rho parameter of particle in [T*m]"}, { "charge", Bunch_charge ,METH_VARARGS,"Set charge(value) or get charge() the charge of particle in e-charge"}, { "macroSize", Bunch_macroSize ,METH_VARARGS,"Set macroSize(value) or get macroSize() the charge of particle in e-charge"}, { "initBunchAttr", Bunch_initBunchAttr ,METH_VARARGS,"Reads and initilizes bunch attributes from a bunch file"}, diff --git a/src/teapot/teapotbase.cc b/src/teapot/teapotbase.cc index cffa337a..0c2b354a 100644 --- a/src/teapot/teapotbase.cc +++ b/src/teapot/teapotbase.cc @@ -231,32 +231,28 @@ void wrapbunch(Bunch* bunch, double length) void kick(Bunch* bunch, double kx, double ky, double kE, int useCharge) { - double charge = +1.0; - if(useCharge == 1) charge = bunch->getCharge(); - double kxc = kx * charge; - double kyc = ky * charge; - double kEc = kE * charge; + //coordinate array [part. index][x,xp,y,yp,z,dE] double** arr = bunch->coordArr(); - if(kxc != 0.) + if(kx != 0.) { for(int i = 0; i < bunch->getSize(); i++) { - arr[i][1] += kxc; + arr[i][1] += kx; } } - if(kyc != 0.) + if(ky != 0.) { for(int i = 0; i < bunch->getSize(); i++) { - arr[i][3] += kyc; + arr[i][3] += ky; } } - if(kEc != 0.) + if(kE != 0.) { for(int i = 0; i < bunch->getSize(); i++) { - arr[i][5] += kEc; + arr[i][5] += kE; } } } @@ -273,6 +269,7 @@ void kick(Bunch* bunch, double kx, double ky, double kE, int useCharge) // i = particle index // pole = multipole number // kl = integrated strength of the kick [m^(-pole)] +// kl already has information about the charge of particle. // skew = 0 - normal, 1 - skew // // RETURNS @@ -282,9 +279,11 @@ void kick(Bunch* bunch, double kx, double ky, double kE, int useCharge) void multpi(Bunch* bunch, int i, int pole, double kl, int skew, int useCharge) { - double charge = +1.0; - if(useCharge == 1) charge = bunch->getCharge(); - double klc = kl * charge; + if(bunch->getCharge() == 0.){ + return; + } + + double klc = kl; std::complex z, zn; double kl1; @@ -326,6 +325,7 @@ void multpi(Bunch* bunch, int i, int pole, double kl, int skew, int useCharge) // pole = multipole number // pole = 0 for dipole, pole = 1 for quad, pole = 2 for sextupole, pole = 3 for octupole // kl = integrated strength of the kick [m^(-pole)] +// kl already has information about the charge of particle. // skew = 0 - normal, 1 - skew // // RETURNS @@ -335,9 +335,11 @@ void multpi(Bunch* bunch, int i, int pole, double kl, int skew, int useCharge) void multp(Bunch* bunch, int pole, double kl, int skew, int useCharge) { - double charge = +1.0; - if(useCharge == 1) charge = bunch->getCharge(); - double klc = kl * charge; + if(bunch->getCharge() == 0.){ + return; + } + + double klc = kl; std::complex z, zn; double kl1; @@ -382,6 +384,7 @@ void multp(Bunch* bunch, int pole, double kl, int skew, int useCharge) // bunch = reference to the macro-particle bunch // pole = multipole number // kl = multipole strength +// kl already has information about the charge of particle. // skew = multipole skew // // RETURNS @@ -391,9 +394,11 @@ void multp(Bunch* bunch, int pole, double kl, int skew, int useCharge) void multpfringeIN(Bunch* bunch, int pole, double kl, int skew, int useCharge) { - double charge = +1.0; - if(useCharge == 1) charge = bunch->getCharge(); - double klc = kl * charge; + if(bunch->getCharge() == 0.){ + return; + } + + double klc = kl; std::complex rootm1 = std::complex(0.0, 1.0); SyncPart* syncPart = bunch->getSyncPart(); @@ -492,6 +497,7 @@ void multpfringeIN(Bunch* bunch, int pole, double kl, int skew, int useCharge) // bunch = reference to the macro-particle bunch // pole = multipole number // kl = multipole strength +// kl already has information about the charge of particle. // skew = multipole skew // // RETURNS @@ -501,9 +507,11 @@ void multpfringeIN(Bunch* bunch, int pole, double kl, int skew, int useCharge) void multpfringeOUT(Bunch* bunch, int pole, double kl, int skew, int useCharge) { - double charge = +1.0; - if(useCharge == 1) charge = bunch->getCharge(); - double klc = kl * charge; + if(bunch->getCharge() == 0.){ + return; + } + + double klc = kl; std::complex rootm1 = std::complex(0.0, 1.0); SyncPart* syncPart = bunch->getSyncPart(); @@ -602,6 +610,7 @@ void multpfringeOUT(Bunch* bunch, int pole, double kl, int skew, int useCharge) // bunch = reference to the macro-particle bunch // length = length of transport // kq = quadrupole field strength [m^(-2)] +// kq already has information about the charge of particle. // // RETURNS // Nothing @@ -610,14 +619,14 @@ void multpfringeOUT(Bunch* bunch, int pole, double kl, int skew, int useCharge) void quad1(Bunch* bunch, double length, double kq, int useCharge) { - double charge = +1.0; - if(useCharge == 1) charge = bunch->getCharge(); - double kqc = kq * charge; - if(kqc == 0.) + if(kq == 0. || bunch->getCharge() == 0.) { drift(bunch,length); return; } + + double kqc = kq; + double x_init, xp_init, y_init, yp_init; double sqrt_kq, kqlength; double cx, sx, cy, sy; @@ -765,6 +774,7 @@ void quad3(Bunch* bunch, double length, double kq, int useCharge) // PARAMETERS // bunch = reference to the macro-particle bunch // kq = strength of quad +// kq already has information about the charge of particle. // // RETURNS // Nothing @@ -773,9 +783,11 @@ void quad3(Bunch* bunch, double length, double kq, int useCharge) void quadfringeIN(Bunch* bunch, double kq, int useCharge) { - double charge = +1.0; - if(useCharge == 1) charge = bunch->getCharge(); - double kqc = kq * charge; + if(bunch->getCharge() == 0.){ + return; + } + + double kqc = kq; double KNL, x_init, xp_init, y_init, yp_init, detM; SyncPart* syncPart = bunch->getSyncPart(); @@ -832,6 +844,7 @@ void quadfringeIN(Bunch* bunch, double kq, int useCharge) // PARAMETERS // bunch = reference to the macro-particle bunch // kq = strength of quad +// kq already has information about the charge of particle. // // RETURNS // Nothing @@ -840,9 +853,11 @@ void quadfringeIN(Bunch* bunch, double kq, int useCharge) void quadfringeOUT(Bunch* bunch, double kq, int useCharge) { - double charge = +1.0; - if(useCharge == 1) charge = bunch->getCharge(); - double kqc = kq * charge; + if(bunch->getCharge() == 0.){ + return; + } + + double kqc = kq; double KNL, x_init, xp_init, y_init, yp_init, detM; SyncPart* syncPart = bunch->getSyncPart(); @@ -1352,14 +1367,12 @@ void bendfringeOUT(Bunch* bunch, double rho) void soln(Bunch* bunch, double length, double B, int useCharge) { //if solenoid field in [T] is zero we have just a drift - if(abs(B) < 1.0e-100){ + if(abs(B) < 1.0e-100 || bunch->getCharge() == 0.){ drift(bunch,length); return; } - - double charge = +1.0; - if(useCharge == 1) charge = bunch->getCharge(); - double Bc = B * charge; + + double Bc = B * bunch->getCharge(); double KNL, phase, cs, sn; double cu, cpu, u_init, pu_init, u, pu, phifac; @@ -1526,6 +1539,8 @@ void RingRF(Bunch* bunch, double ring_length, int harmonic_numb, } SyncPart* syncPart = bunch->getSyncPart(); + + double p_synch_in = syncPart->getMomentum(); if(phase_s != 0.) { @@ -1533,6 +1548,10 @@ void RingRF(Bunch* bunch, double ring_length, int harmonic_numb, kin_e += coeff * voltage * sin(phase_s); syncPart->setMomentum(syncPart->energyToMomentum(kin_e)); } + + double p_synch_out = syncPart->getMomentum(); + + double xp_yp_coeff = p_synch_in/p_synch_out; //coordinate array [part. index][x,xp,y,yp,z,dE] double** arr = bunch->coordArr(); @@ -1541,6 +1560,9 @@ void RingRF(Bunch* bunch, double ring_length, int harmonic_numb, { deltaV = voltage * ( sin(harmonic_numb*Factor*arr[i][4] + phase_s)); arr[i][5] += coeff * deltaV; + + arr[i][1] *= xp_yp_coeff; + arr[i][3] *= xp_yp_coeff; } } diff --git a/src/teapot/wrap_teapotbase.cc b/src/teapot/wrap_teapotbase.cc index fe8a29ca..b353130d 100644 --- a/src/teapot/wrap_teapotbase.cc +++ b/src/teapot/wrap_teapotbase.cc @@ -488,7 +488,7 @@ extern "C" {"rotatexy", wrap_rotatexy, METH_VARARGS, "Rotates bunch around z axis "}, {"drift", wrap_drift, METH_VARARGS, "Tracking a bunch through a drift "}, {"drifti", wrap_drifti, METH_VARARGS, "Drifts one macroparticle in the bunch"}, - {"wrapbunch", wrap_wrapbunch, METH_VARARGS, "Tracking a bunch through a wrapbunch routine"}, + {"wrapbunch", wrap_wrapbunch, METH_VARARGS, "Tracking a bunch through a wrapbunch routine"}, {"multp", wrap_multp, METH_VARARGS, "Tracking a bunch through a multipole "}, {"multpfringeIN", wrap_multpfringeIN, METH_VARARGS, "Tracking a bunch through an IN edge of a multipole "}, {"multpfringeOUT", wrap_multpfringeOUT, METH_VARARGS, "Tracking a bunch through an OUT edge of a multipole"},