ocelot-collab / ocelot

OCELOT is a multiphysics simulation toolkit designed for studying FEL and storage ring-based light sources.
GNU General Public License v3.0
85 stars 58 forks source link

Bend types #242

Closed koesterf closed 5 months ago

koesterf commented 5 months ago

Hello friends,

I am looking for a standard "RBend" C-type chicane example but couldn't found anything. Is there any example for RBend ?

Cheers, Frank

st-walker commented 5 months ago

Dear Friend,

It exists:

from ocelot.cpbd.elements import RBend
rb = RBend(l=1, angle=0.4) # etc.

Best wishes, Stuart

koesterf commented 5 months ago

Dear Stuart,

Thanks for being so fast ! I am trying to understand how RBend, SBend and Bend are implemented in the code. So far I haven't seen any caution for the user about magnet length. For instance, for SBend it is the effective length that one needs to take care (length that is modified by the bending angle) etc.

Is everything is taken care inside the code structure or the user have to take care of these things ?

It might be a good idea to include some examples of all kinds of bending magnets and their ability in beam transportation and manipulation.

Cheers, Frank

st-walker commented 5 months ago

Yeah the magnets are a bit odd in their implementation and "it" (i.e. the implementation) is always defined somewhere else.

I haven't thought about it but I believe it must follow the mad8 convention. So that would mean the length is the chord length. So it has no meaning to "take care" of the length/ange combination in this case, as they are independent of each other. @sergey-tomin is this correct?

koesterf commented 5 months ago

@st-walker,

If I look at Tutorial N1 and Tutorial N7, then I could see that "Bend" type magnet is used. In tutorial N1, the bend is defined as

# defining of the bending magnet
B = Bend(l=2.7, k1=-.06, angle=2*pi/16., e1=pi/16., e2=pi/16.)

including k1 and term angle is different from entrance (e1) and exit (e2) angle.

Whereas, in tutorial 7, the bend is defined very similar to SBend

b1 = Bend(l = 0.5, angle=-0.0336, e1=0.0, e2=-0.0336, gap=0, tilt=0, eid='BB.393.B2')
b2 = Bend(l = 0.5, angle=0.0336, e1=0.0336, e2=0.0, gap=0, tilt=0,  eid='BB.402.B2')
b3 = Bend(l = 0.5, angle=0.0336, e1=0.0, e2=0.0336, gap=0, tilt=0, eid='BB.404.B2')
b4 = Bend(l = 0.5, angle=-0.0336, e1=-0.0336, e2=0.0, gap=0, tilt=0, eid='BB.413.B2')

but now with effective drift length

d = Drift(l=1.5/np.cos(b2.angle))

Effective length consideration appears more accurate, however, such kind of things makes the user confused.

How to decide what the user should do ?

The user should consider only the physical values of the magnet or the user should control everything carefully by themselves.

koesterf commented 5 months ago

Dear @st-walker, @sergey-tomin

I have a question for both of you. I made a schematic of RBend magnetic chicane. Please, see the attached pdf file. If you could see and answer my question raised there then I may get my answer about RBend.

RBend.pdf

st-walker commented 5 months ago

No that's wrong, or at least not how we do it in euxfel: e.g.

bl_48i_i1 = SBend(l=0.2, angle=-0.099484, e2=-0.099484, eid='BL.48I.I1')
bl_48ii_i1 = SBend(l=0.2, angle=0.099484, e1=0.099484, eid='BL.48II.I1')
bl_50i_i1 = SBend(l=0.2, angle=0.099484, e2=0.099484, eid='BL.50I.I1')
bl_50ii_i1 = SBend(l=0.2, angle=-0.099484, e1=-0.099484, eid='BL.50II.I1')

note alternating use of e2 and e1.

idk anything about Bend, tbh i thought it was just a base class, not to be instantiated. i only ever use sbend personally.

drifts are the true path length of the bunch, NOT the projected length. so should be l1/cos(theta) NOT l1 (referring 2 your diagram). so yes u got that right in ur diagram. in summary i think ur diagram is correct minus any considerations re e1/e2.

koesterf commented 5 months ago

@st-walker,

in EUXFEL beamline simulations do you use SBend for rectangular magnets ?

For the case shown in my diagram (rectangular magnets), can I use SBend ? OR rectangular magnets should only be simulated using RBend.

st-walker commented 5 months ago

rbend is not what you want. the incident angle to RBEND is half the bend angle, in a chicane the incident angle is (typically) 0 degrees.. if you want all 4 faces of dipoles to be parallel (like in your diagram), you need sbend with e1 and e2 set like in the piece of code in my previous comment. to be clear, an rbend is an sbend with appropriately set e1 and e2. they are not so different.

in EUXFEL beamline simulations do you use SBend for rectangular magnets ?

yes.

For the case shown in my diagram (rectangular magnets), can I use SBend ?

yes, but set e1 and e2 correctly.

koesterf commented 5 months ago

Dear @st-walker,

Thank you for correcting my mis-conception.

Can you help me to figure out if a 4 rectangular magnets based chicane (as shown in my figure) can be simulated with the SBend as given here:

dc_12 = Drift(1.0)  # distance between chicane 1 & 2     (c is for chicane)
dc_23 = Drift(0.6)  # distance between chicane 2 & 3
dc_34 = Drift(1.0)  # distance between chicane 3 & 4

Lmag = 0.80         # length of the magnets

angle0 = theta

# converting to actual path
dc_12 = dc_12 / np.cos(theta) 
dc_34 = dc_34 / np.cos(theta) 

# magnets configuration
b1 = SBend(l=Lmag, angle=theta, e1=0.0, e2=theta, eid='B1')
b2 = SBend(l=Lmag, angle=-theta, e1=-theta, e2=0.0, eid='B2')
b3 = SBend(l=Lmag, angle=-angle0, e1=0.0, e2=-theta, eid='B3')
b4 = SBend(l=Lmag, angle=angle0, e1=theta, e2=0.0, eid='B4')

chicane = (b1, dc_12_actual , b2, dc_23, b3, dc_34, b4)

I assume that for each rectangular magnet e1 is the entrance angle and e2 is the exit angle.

I see the C-type rectangular magnet based chicane as follows: entrance angle of 1st chicane is zero and exit angle is theta, similarly the entrance angle of 2nd chicane is -(theta) and exit angle is zero, the entrance angle of 3rd chicane is zero and exit angle is -(theta), the entrance angle of 4th chicane is theta and exit angle is zero.

Is it correct ?

st-walker commented 5 months ago

I assume that for each rectangular magnet e1 is the entrance angle and e2 is the exit angle

they are the poleface rotations for the entrance (e1) and exit(e2). for an sbend, 0 means beam enters (exits) perpendicular to the magnet face.

Is it correct ?

looks good to me i think, this way all 8 faces are parallel, like in your diagram. if you did it like in https://github.com/ocelot-collab/ocelot/issues/242#issuecomment-2096895341 then it is correct.

except:

># converting to actual path
dc_12 = dc_12 / np.cos(theta) 
dc_34 = dc_34 / np.cos(theta) 

Looks fine, but you don't actually use these values in your snippet you shared. You are still using the projected lengths, whereas you need the path lengths for the drifts. indeed i did not run this but it looks like it will crash because you are dividing a Drift by a float.

koesterf commented 5 months ago

Dear @st-walker,

Thanks for making things even more clear. OK, so as far as I understood, you think that the way I am using SBend for rectangular magnets based chicane is correct except the length conversion part dc_12 = dc_12 / np.cos(theta).

So, I should simply do as:

L1 = 1.0
L2 = 0.6
L3 = 1.0

dc_12 = Drift(L1)   # distance between chicane 1 & 2     (c is for chicane)
dc_23 = Drift(L2)   # distance between chicane 2 & 3
dc_34 = Drift(L3)   # distance between chicane 3 & 4

Lmag = 0.80         # length of the magnets

angle0 = theta

# converting to actual path
#dc_12_actual = Drift( L1 / np.cos(theta)  )
#dc_34_actual = Drift( L3 / np.cos(theta) )

# magnets configuration
b1 = SBend(l=Lmag, angle=theta, e1=0.0, e2=theta, eid='B1')
b2 = SBend(l=Lmag, angle=-theta, e1=-theta, e2=0.0, eid='B2')
b3 = SBend(l=Lmag, angle=-angle0, e1=0.0, e2=-theta, eid='B3')
b4 = SBend(l=Lmag, angle=angle0, e1=theta, e2=0.0, eid='B4')

#chicane = (b1, dc_12_actual , b2, dc_23, b3, dc_34_actual, b4)
chicane = (b1, dc_12 , b2, dc_23, b3, dc_34, b4)

which chicane part is correct ?

chicane = (b1, dc_12_actual , b2, dc_23, b3, dc_34_actual, b4) OR chicane = (b1, dc_12 , b2, dc_23, b3, dc_34, b4)

st-walker commented 5 months ago

looks reasonable this one is correct, as i said:

chicane = (b1, dc_12_actual , b2, dc_23, b3, dc_34_actual, b4)

koesterf commented 5 months ago

Dear @st-walker,

Thank you for clearing my mis-conception. You indeed are great !

Cheers, Frank

koesterf commented 5 months ago

Dear @st-walker

what is the unit of chirp in the code. Normally, chirp is defined as per-meter. In ocelot the chirp is normalized ? Because when I use chirp = 0.01 it gives me quite a large band width.

If chirp = 0.01 what is its unit ?

Cheers, Frank

st-walker commented 5 months ago

which function are you talking about? I can only guess

in generate_parray, please read the docstring:

:param chirp: energy chirp [unitless], linear correlation - p_i += chirp * tau_i/sigma_tau

koesterf commented 5 months ago

@st-walker, I want to calculate (1 + h*R56), here h is the energy chirp of the electron beam and R56 is the longitudinal dispersion. R56 can be calculated directly but without knowing the chirp (h) in practical units, it is difficult to calculate (1 + h*R56) correctly.

st-walker commented 5 months ago

where are you getting h from in ocelot though? Where is this value of chirp that you have in the wrong units from ocelot?

koesterf commented 5 months ago

@st-walker this is what I was wondering about .... how to get the value of h from ocelot !

In literature, the equation 1 + h*R56 is used as a baseline for obtaining maximum compression in chicane. For good chicane design its value should be negative but as close to zero as possible. That is the reason I want to check 1 + h*R56. Any help will be very appreciable.

st-walker commented 5 months ago

this is what I was wondering about .... how to get the value of h from ocelot

OK but that's not what you asked originally, that is a new question.

How would I do it? ocelot.cpbd.beam.global_slice_analysis(parray) to get the mean slice energies (among aother things), and then fit a function to the slice s vs energy and then you get all the chirps you want. I've never calculated the chip from a beam in OCELOT. But I'd try that.

sergey-tomin commented 5 months ago

Hi As far as I understood, you want to generate a beam with a certain chirp, let's say h = 0.1 GeV/m. @st-walker already told you that the generate_parray() function uses the following expression for chirp: p_i += chirp * tau_i/sigma_tau, and I think it is enough to figure out how to do it.

I would write chirp = h / E_beam (in GeV) * sigma_tau, please double-check it. This is not a precise equation since in Ocelot, p is dE/pc, but it is a good approximation for relativistic beams.

koesterf commented 5 months ago

Thanks @st-walker and @sergey-tomin

I will follow your suggestions.

By the way, I also thought how to do it and come up with a simple way (but not sure if this is correct or not). As h=(1/E)*(dE/dz), E is the mean energy of the beam, dE is the energy spread and dz corresponding beam length. Maybe, using a longitudinal phase-space, I can calculate the FWHM of energy spread and calculate the corresponding dz and then calculate h. But first, I will try your methods.

Cheers, Frank

st-walker commented 5 months ago

You can do it that way. Your way is easier but it would include uncorrelated energy spread in chirp calculation. For very large chirp compared to uncorrelated energy spread this is fine I guess. I think you should try your way first since it is easy., just parray.p().std() / parray.tau().std() = h. p is already normalised w.r.t E (or p, does not matter for high energy). so this is fine.

koesterf commented 5 months ago

Thanks @st-walker I guess now I can do what I was looking for. Thanks for your great help.