GitLab will be upgraded on June 23rd evening. During the upgrade the service will be unavailable, sorry for the inconvenience.

Commit 5c0c50dd authored by Emmanuel Papillon's avatar Emmanuel Papillon

add undulator tool to comptue undulator position from energy

parent 50e04745
import numpy as np, pylab as pl
# Code to fit / predict undulator peak positions
# A script by Jon Wright, 2018
#
# 0.95 Er[GeV]^2
# Ef[keV] = ---------------------
# ( 1 + k^2 / 2 ) p[cm]
#
# Emax[keV] is the maximum energy for the fundamental at k=0
#
# Emax[keV] = 0.95 Er[GeV]^2 / p[cm]
#
#
# Using this we have
#
# Emax[keV]
# Ef[keV] = ---------------
# ( 1 + k^2 / 2 )
#
#
# 0.95 Er[GeV]^2
# k^2/2 = --------------- - 1 = Emax/Ef - 1
# Ef[keV] p[cm]
#
#
#
# k is the deflection parameter
# k ~ field ~ exp(-gap)
#
# log(k) vs gap should be more-or-less linear
# log(k^2) = 2 log(k)
# log(k^2/2) = 2 log(k) - 2 log(2)
# ... so we will work with Emax and k^2/2 in the fitting
class undulator( object ):
def __init__(self, period=1.8, Er=6., pars=[0,0],
mingap = 5.9, maxgap = 30):
"""
period is in cm
Er is the storage ring energy [GeV]
pars = parameters for the model used here (gap vs log(Emax/Ef - 1))
polynomial coefficients (order = 1)
mingap = minimum allowed gap [mm]
maxgap = maximum allowed gap [mm]
"""
self.period = period
self.pars = pars
self.Er = Er
self.Emax = 0.95*self.Er*self.Er/self.period
self.mingap = mingap
self.maxgap = maxgap
def __repr__(self):
""" Something to copy/paste to create one """
s = "undulator( period = %f, Er = %f, pars = [%f, %f] )"%(
self.period, self.Er, self.pars[0], self.pars[1] )
return s
def k_from_gap(self, gap):
"""
Compute k for a given gap [mm]
"""
assert self.pars is not None
log_kk = np.polyval( self.pars, np.asarray(gap) )
k = np.sqrt( 2*np.exp( log_kk ) )
return k
def ef_from_gap(self, gap):
"""
Compute the energy [keV] of the fundamental for a given gap [mm]
"""
assert self.pars is not None
kk = np.exp( np.polyval( self.pars, np.asarray(gap) ) )
return self.Emax / (1 + kk)
def gap_from_ef(self, ef):
"""
Compute the gap [mm] needed for get a specific energy [keV]
"""
assert ef <= self.Emax
log_kk = np.log( self.Emax/np.asarray(ef) - 1 )
# log_kk = p[0]*g + p[1]
g = (log_kk - self.pars[1])/self.pars[0]
return g
def k_from_ef( self, ef ):
"""
Compute the k given the fundamental energy [keV]
"""
kk = self.Emax / np.asarray(ef) - 1
assert (kk >= 0).all()
k = np.sqrt(2*kk)
return k
def fitgaps( self, gaps, ef ):
"""
Fit undulator peak position from a gap scan (or scans)
g = fitted gap values [mm]
ef = energy of fundamental [keV] (e.g. monochromator energy/harmonic)
Uses Emax = 0.95 * Er[GeV]^2 / period[cm]
"""
kk = self.Emax/np.asarray(ef) - 1 # k^2/2
self.pars = np.polyfit( gaps, np.log(kk), 1 ) # fits a line
def plotfit( self, gaps, ef):
"""
Diagnostic plots to see if a fit is OK
"""
gc = np.linspace( self.mingap, self.maxgap, 100)
kcalc = self.k_from_gap( gc )
kobs = self.k_from_ef( ef )
pl.figure()
pl.subplot(121)
pl.plot( gaps, kobs, "o")
pl.plot( gc, kcalc, "-")
pl.semilogy()
pl.ylabel( "k")
pl.xlabel("gap, mm")
pl.title("period %.2f"%(self.period))
pl.subplot(122)
pl.plot( gaps, ef, "o")
efcalc = self.ef_from_gap( gc )
pl.plot( gc, efcalc, "-")
pl.xlabel("gap, mm")
pl.ylabel("Ef[keV]")
def peaks( self, energy ):
"""
Return a list of peaks (n, gap) for
given energy
"""
i = 0
hg = []
while 1:
i += 1
ef = energy / i
if ef < self.Emax:
g = self.gap_from_ef( ef )
if g < self.mingap:
break
hg.append( (i,g) )
return hg
def test_u22():
# u22:
e = 65.351
n = 9., 7., 5.
g = 7.051, 8.709, 12.80
ef = [ e/ni for ni in n ]
u22 = undulator( period = 2.2 )
u22.fitgaps( g, ef )
u22.plotfit( g, ef )
print( "# u22 = ", u22 )
for e in [40, 65.351, 78, 88]:
print(e," keV:")
print( u22.peaks( e ) )
def test_cpm18():
# cpm18
e = 65.351
n = 7., 5.
g = 6.748, 8.931
ef = [ e/ni for ni in n ]
cpm18 = undulator( period = 1.8 )
cpm18.fitgaps( g, ef )
cpm18.plotfit( g, ef)
print( "# cpm18 = ", cpm18 )
for e in [40, 65.351, 78, 88]:
print( e," keV:")
print( cpm18.peaks( e ) )
print("k at 6mm",cpm18.k_from_gap(6))
def ebs_cpm18():
e1 = 65.3508 # Hf
pks = [6.057, 6.758, 7.68, 9.0, 12.06]
nhf = [8,7,6,5,4]
ef = [e1/ni for ni in nhf]
epb=88.0045
pks+= [6.448, 7.028, 7.700, 8.676, 10.28]
npb=[10, 9, 8, 7, 6]
ef +=[epb/ni for ni in npb]
cpm18 = undulator( period = 1.8 )
cpm18.fitgaps( pks, ef )
cpm18.plotfit( pks, ef )
print(cpm18.pars)
def ebs_u22():
e1 = 65.3508 # Hf
pks = 7.111, 7.858, 8.800, 10.32, 13.05
n = 9, 8, 7, 6, 5
u22 = undulator( period = 2.2 )
ef = [e1/ni for ni in n]
u22.fitgaps( pks, ef )
u22.plotfit( pks, ef )
print(repr(u22.pars))
def main():
test_u22()
test_cpm18()
ebs_cpm18()
ebs_u22()
pl.show()
id11_cpm18 = undulator( period = 1.800000, Er = 6.000000,
pars = [-0.35011249 , 2.3931669 ] )
# pars = [-0.377876, 2.584470] )
id11_u22 = undulator( period = 2.200000, Er = 6.000000,
pars = [-0.30146691, 2.26371134])
# pars = [-0.311363, 2.317529] )
def undpeaks( en ):
print("Energy = %.3f keV"%(en))
print(" N gap/mm cpm18")
for i, g in id11_cpm18.peaks(en):
print("%2d %6.3f"%(i,g))
print(" N gap/mm u22")
for i, g in id11_u22.peaks(en):
print("%2d %6.3f"%(i,g))
def __main__():
import sys
undpeaks(float(sys.argv[1]))
__all__=[ undpeaks , ]
#if __name__=="__main__":
# Not allowed for user_scripts...
# print("I think I am in __main__")
#__main__()
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment