The OPAL Discussion Forum

["Snuverink Jochem (PSI)" <jochem.snuverink AT psi.ch>]

Hi Finn,

Many thanks for the detailed email.

Unless I misunderstood something, I think you discovered two bugs in Opal-T.

1) qi_m, the charge per macroparticle is defined in TrackRun.cpp (last line):

    // Return charge per macroparticle.

    return beam->getCharge() * beam->getCurrent() / (beam->getFrequency()*1.0e6) / numberOfParticles;

The factor beam->getCharge() shouldn't be there.

2) As you discovered, the charge (only its sign) is not taken into account when calculating the bending field (in the 1DPROFILE1-DEFAULT case). This is the case both for RBend and SBend.

I'll make two issues in the issue tracker.

As a workaround for each case:

1) You can halve the current, since the current only used to calculate the bunch charge.

2) You can try to reduce the design energy of the bend such to halve beta*gamma:

new_design_energy = m0 / 2 * ( sqrt ( ( design_energy / m0 + 1 ) ^ 2 + 3 ) - 1 )

Best regards,

Jochem

On 22 May 2020, at 02:20, finn.oshea AT nusano.com wrote:

I've been trying to get a simulation going with alpha particles and the
various printouts suggest that charge isn't being handled consistently in
Opal-t.

I had a long version of this message with more output from Opal-t, but sympa
would not let me send it.

Starting with a FODO example that Andreas sent me, I can change the charge of
the bunch by changing the charge of the particles, in this case simply
changing from protons to carbon atoms increases the bunch charge by a factor
of 12, even though the current has not been changed.

Since the total bunch charge is defined through the relationship BCURRENT = Q
* BFREQ and neither BCURRENT nor BFREQ has changed, Q shouldn't either, but it
does.

This problem also seems to have crept in to the computation of the field
strength of an SBEND in a way that seems doubly incorrect to me.

Here is the input to Opal-t for a particle of charge +2 (I've deleted some
lines that do not matter, this is not the same beam as shown above):
OPAL>
OPAL> value: {EDES,GAMMA,BETA,P0,BRHO}
={0.05064,1.01359,0.163181,0.616503,1.02822}
OPAL>

  54 REAL switch_angle = 15.0 * Pi / 180.0; // 15-degree switch angle
  55 REAL switch_radius = 1; // 1 meter bending radius
  56
  57 MDL_S.01 : SBEND, ANGLE = switch_angle,
  58           FMAPFN = "1DPROFILE1-DEFAULT",
  59           ELEMEDGE = 2.664,
  60           DESIGNENERGY = design_energy,
  61           L = 2 * switch_radius * sin(switch_angle / 2.0),
  62           E1 = 0, E2 = 0, // make a "true" sector magnet

P0 and BRHO give the particle momentum and rigidity.

However, when I look at the output of the SBEND magnet (info level 2) I see
that the strength of the magnet is too large by a factor of 2:
SBend [2]> Reference Trajectory Properties
SBend [2]>
======================================================================
SBend [2]>
SBend [2]> Bend angle magnitude:    0.261799 rad (15 degrees)
SBend [2]> Entrance edge angle:     0 rad (0 degrees)
SBend [2]> Exit edge angle:         0 rad (0 degrees)
SBend [2]> Bend design radius:      1 m
SBend [2]> Bend design energy:      5.0655e+07 eV
SBend [2]>
SBend [2]> Bend Field and Rotation Properties
SBend [2]>
======================================================================
SBend [2]>
SBend [2]> Field amplitude:         2.05674 T
SBend [2]> Field index:  0
SBend [2]> Rotation about z axis:   0 rad (0 degrees)
SBend [2]>
SBend [2]> Reference Trajectory Properties Through Bend Magnet with Fringe
Fields
SBend [2]>
======================================================================
SBend [2]>
SBend [2]> Reference particle is bent: 0.2618 rad (15 degrees) in x plane
SBend [2]> Reference particle is bent: 0 rad (0 degrees) in y plane

Where I've deleted some lines that I don't think matter.  Changing the charge
of the particle does not change the output.  The field strength is still
2.05674 T.  When I look in to the code that produces this text (Bend2D.cpp) it
appears to get that value from the following function (line 1070 of the file):

void Bend2D::setBendStrength() {

   // Estimate bend field magnitude.
   double mass = RefPartBunch_m->getM();
   double gamma = designEnergy_m / mass + 1.0;
   double betaGamma = std::sqrt(std::pow(gamma, 2.0) - 1.0);
   double charge = RefPartBunch_m->getQ();

   fieldAmplitude_m = ((charge / std::abs(charge)) * betaGamma * mass /
                       (Physics::c * designRadius_m));


   // Search for angle if initial guess is not good enough.
   findBendStrength(mass, gamma, betaGamma, charge);
}

This function is interesting for two reasons:
(1) The formula being used is q*P0 / (rho*c) where it should be P0/(q*rho*c)
(2) For some reason the charge is divided by its absolute value, perhaps the
author just wanted the sign of the charge
This formula will always return the wrong value for particles of charge
greater than 1 (in an absolute sense).

Nevertheless, I tracked down the getQ() function to the file PartBunchBase.hpp
(which also happens to be the source of the BUNCH output producing the
erroneous bunch charge value), so perhaps that is the root of the problem.
But it does appear that the incorrect charge output is not just a superficial
print bug, but rather the value for the bunch charge really is incorrect in
the variable qi_m.  I can't figure out where qi_m comes from, I've never used
C++.

Thanks for any assistance,
Finn