Adding slope errors to CRLs (double paraboloid lens) about xrt

[Faryzww]

I want to know if there's something wrong of this code, and if it is right for us to add the slope errors to CRLs like this.
Best regards.

[kklmn]

Hello,
Where should I look at? Which part of the code is yours? Where in the code "the slope error is R=100um"? What is "CRL’s R" and where in the code you use it?

[Faryzww]

I upload my codes. CRL_ERR_TEST_Beamline.py is the main code, and Distorted_Element.py contains the function of distorted double paraboloid lens.
The initial radius of CRLs is set to be 0.1mm. In order to know whether my code works right, I take a CRL's (R=0.1mm) slope as an error and add it to the initial CRLs. The result focal length should be that of CRLs with R=0.05mm. But I failed. Moreover, the spot after adding slope error is very strange with a dark in the middle. I am so confused.

[Faryzww]

Distorted_Element.py contains:

-- coding: utf-8 --

"""
Created on Thu Jul 26 11:55:05 2018
2018.8.15:修改

"""

import os, sys;

sys.path.append(os.path.join('..', '..', '..')) # analysis:ignore
import numpy as np
import matplotlib.pyplot as plt

import re

from scipy import ndimage

from scipy import optimize

import xrt.backends.raycing as raycing
import xrt.backends.raycing.oes as roe
import xrt.backends.raycing.materials as rm

===================================================================================================#

class DoubleParaboloidLens_Distorted(roe.DoubleParaboloidLens):

所有的self代表的是class的对象

hiddenMethods = roe.DoubleParaboloidLens.hiddenMethods + ['double_refract']  # ?

def __init__(self, *args, **kwargs):
    kwargs = self.__pop_kwargs(**kwargs)
    roe.DoubleParaboloidLens.__init__(self, *args, **kwargs)
    self.read_data()

def __pop_kwargs(self, **kwargs):
    self.fname1 = kwargs.pop('fname1')  # fname1 and fname2 is the ANSYS data file, the mirror surface profile.
    self.fname2 = kwargs.pop('fname2')
    self.fname3 = kwargs.pop('fname3')
    self.get_distorted_surface = kwargs.pop('get_distorted_surface')
    return kwargs

def assign_auto_material_kind(self, material):
    material.kind = 'lens'

def read_data(self):
    if self.get_distorted_surface is None:
        return
    self.warpX, self.warpY, self.warpZ = Read_distorted_surface(self.fname1, self.fname2, self.fname3)
    self.warpNX, self.warpNY = len(self.warpX), len(self.warpY)
    self.limPhysX0 = np.min(self.warpX), np.max(self.warpX)
    self.limPhysY0 = np.min(self.warpY), np.max(self.warpY)
    self.get_surface_limits()
    self.warpA, self.warpB = np.gradient(self.warpZ)
    dx = self.warpX[1] - self.warpX[0]
    dy = self.warpY[1] - self.warpY[0]
    self.warpA = np.arctan(self.warpA / dx)
    self.warpB = np.arctan(self.warpB / dy)

    self.warpSplineZ = ndimage.spline_filter(self.warpZ)
    self.warpSplineA = ndimage.spline_filter(self.warpA)
    self.warpSplineB = ndimage.spline_filter(self.warpB)

def local_z1(self, x, y):
    """Determines the normal vector of OE at (x, y) position."""
    z = (x ** 2 + y ** 2) / (4 * self.focus)
    coords = np.array(
        [(x - self.limPhysX0[0]) /
         (self.limPhysX0[1] - self.limPhysX0[0]) * (self.warpNX - 1),
         (y - self.limPhysY[0]) /
         (self.limPhysY0[1] - self.limPhysY0[0]) * (self.warpNY - 1)])

    z = z + ndimage.map_coordinates(self.warpSplineZ, coords, prefilter=True)
    if self.zmax is not None:
        z[z > self.zmax] = self.zmax
    return z

def local_z2(self, x, y):
    """Determines the surface of OE at (x, y) position."""

认为第二个面和第一个面一样

    return self.local_z1(x, y)


def local_n1(self, x, y):
    Ax = -self.warpSplineA
    By = -self.warpSplineB
    coords = np.array(
        [(x - self.limPhysX0[0]) /
         (self.limPhysX0[1] - self.limPhysX0[0]) * (self.warpNX - 1),
         (y - self.limPhysY0[0]) /
         (self.limPhysY0[1] - self.limPhysY0[0]) * (self.warpNY - 1)])
    a0 = ndimage.map_coordinates(Ax, coords, prefilter=True)
    b0 = ndimage.map_coordinates(By, coords, prefilter=True)
    if self.zmax is not None:
        z = (x ** 2 + y ** 2) / (4 * self.focus)
        z = z + ndimage.map_coordinates(self.warpSplineZ, coords, prefilter=True)

        a = -x / (2 * self.focus) + a0  # -dz/dx  含误差部分
        b = -y / (2 * self.focus) + b0  # -dz/dy  含误差部分

        if isinstance(a, np.ndarray):
            a[z > self.zmax] = 0
        if isinstance(b, np.ndarray):
            b[z > self.zmax] = 0

        c = np.ones_like(x)
        norm = (a ** 2 + b ** 2 + 1) ** 0.5
    return [a / norm, b / norm, c / norm]

def local_n2(self, x, y):
    return self.local_n1(x, y)

if name == "main":
see_the_bump()

[Faryzww]

And my main code of CRL is as follows(CRL_ERR_TEST_Beamline.py):

import os, sys; sys.path.append(os.path.join('..', '..', '..')) # analysis:ignore
import os, sys; sys.path.append(os.path.join('C:\Users\Lenovo\Desktop\XRT_Simulation\surface error\Error_Generation'))

import numpy as np
import matplotlib.pyplot as plt
#import matplotlib as mpl
#mpl.use('agg')

import xrt.backends.raycing as raycing
import xrt.backends.raycing.sources as rs
import xrt.backends.raycing.oes as roe
import xrt.backends.raycing.apertures as ra
import xrt.backends.raycing.run as rr
import xrt.backends.raycing.materials as rm
import xrt.backends.raycing.screens as rsc
import xrt.plotter as xrtp
import xrt.runner as xrtr

import Distorted_Element as Distorted

showIn3D = False

Lens = roe.DoubleParaboloidLens
Lens_Err=Distorted.DoubleParaboloidLens_Distorted

R = 0.1
focus = R/2

mm，R=2*focus

ztatal = 0.255
t=0.03
zmax = (ztatal - t)/2

dz = 3.

E0 = 12400.

eV

f = 1770.

mm

p = 50000
print('q = {0:.2f}'.format(1/(1/f-1/p)))

mm

xyLimits = -5, 5

mBeryllium = rm.Material('Be', rho=1.845, kind='lens')
material = mBeryllium

crystalSi01 = rm.CrystalSi(name=None)

dE = E0*0.0001/2
eMin0,eMax0 = E0-dE,E0+dE

nrays = 1e4
kwargs_B4 = dict(
nrays=nrays, center=[0, 0, 0],
eE=6.0, eI=0.2, eEspread=0.00111,
eEpsilonX=0.02728, eEpsilonZ=0.00276,
betaX=10.12, betaZ=9.64,
xPrimeMax=0.0012, zPrimeMax=0.0012,

xPrimeMaxAutoReduce=True, zPrimeMaxAutoReduce=True,

    targetE=[E0, 3], eMin=eMin0, eMax=eMax0, 
    uniformRayDensity=True,#Faulse        
    K=r"auto",period=32,n=156) 

#p_slitWM = 29000
#anglex_WM,anglez_WM = 10e-6,10e-6#27.79e-6,25.33e-6
#slitDxWM,slitDzWM = anglex_WMp_slitWM, anglez_WMp_slitWM

def build_beamline(nrays=1e4):

beamLine = raycing.BeamLine()    

beamLine.source = rs.Undulator(beamLine,**kwargs_B4)

beamLine.fsm1 = rsc.Screen(beamLine, 'FSM1', (0, p - 500, 0))



beamLine.roundAperture01 = ra.RoundAperture(
    bl=beamLine,
    name=None,
    center=[0, 49999, 0],
    r=0.06)


beamLine.lenses = []     
lens_error = 1#0-no error，1-with error

Surface_name='CRL_try_'   
kwargs_OE1 = dict(pitch=np.pi/2, t=t, material=material,
                limPhysX=[-0.12/2, 0.12/2], limPhysY=[-0.12/2, 0.12/2], shape='round',
                focus=focus, zmax=zmax, alarmLevel=0.1)      
ilens = 0
while True: 
    center=[0, p + dz*ilens, 0]           
    if lens_error==0:
        lens = Lens(bl=beamLine,name = 'Lens{0:02d}'.format(ilens), 
                center=center,**kwargs_OE1)           
    else:
        kwargs_OE1['get_distorted_surface'] ='error'
        kwargs_OE1['fname1'] =Surface_name+'X.txt'
        kwargs_OE1['fname2'] =Surface_name+'Y.txt'
        kwargs_OE1['fname3'] =Surface_name+'Z.txt'              
        lens = Lens_Err(bl=beamLine,name ='Lens{0:02d}'.format(ilens), 
                center=center,**kwargs_OE1)                  
    beamLine.lenses.append(lens)
    if ilens == 0:
        nCRL = lens.get_nCRL(f, E0)           
    ilens += 1
    if nCRL - ilens < 0.5:  
        break 
    
beamLine.fsm2 = rsc.Screen(beamLine,'FSM2')
beamLine.fsm2.dqs = np.linspace(-1800, 100, 11)  

return beamLine

def run_process(beamLine):

waveOnSlit = beamLine.slitWM.prepare_wave(beamLine.source,nrays)

beamSource = beamLine.source.shine(wave=waveOnSlit,withAmplitudes=True,

fixedEnergy=False)

beamSource = beamLine.source.shine()

outDict = {'beamSource': beamSource} 

beamFSM1 = beamLine.fsm1.expose(beamSource)
outDict['beamFSM1'] = beamFSM1

beamIn = beamSource
for ilens, lens in enumerate(beamLine.lenses):
    lglobal, llocal1, llocal2 = lens.double_refract(beamIn, needLocal=True)
    beamIn = lglobal
    strl = '_{0:02d}'.format(ilens)
    outDict['beamLensGlobal'+strl] = lglobal
    outDict['beamLensLocal1'+strl] = llocal1
    outDict['beamLensLocal2'+strl] = llocal2
    
if Lens==roe.DoubleParaboloidLens:
    nFactor = 0.5
elif Lens==roe.ParabolicCylinderFlatLens:
    nFactor = 2.
else:
    nFactor = 1.
delta=(1. - material.get_refractive_index(E0).real)     
f_real= 2 * focus / len(beamLine.lenses) /delta*nFactor    
q=1/(1/f_real-1/p)
print('Focal length of the lenses ={0:3} mm'.format(f_real))
print('The number of the lenses ={0:3} '.format(ilens))
print('The number of the lenses ={0:3} '.format(ilens))


for i, dq in enumerate(beamLine.fsm2.dqs):        
    beamLine.fsm2.center[1] = p + q + dq
    outDict['beamFSM2_{0:02d}'.format(i)] = beamLine.fsm2.expose(lglobal)

if showIn3D:
    beamLine.prepare_flow() 
          
return outDict

rr.run_process = run_process

def define_plots(beamLine):
plots = []
xrtp.yTextPosNraysR = 0.82
xrtp.yTextPosNrays1 = 0.52
xyLimits = 50
plot00 = xrtp.XYCPlot(
beam=r"beamFSM1",
xaxis=xrtp.XYCAxis(
label=r"x",unit=r"$\mu$m",limits=[-xyLimits,xyLimits]),
yaxis=xrtp.XYCAxis(
label=r"z",unit=r"$\mu$m",limits=[-xyLimits,xyLimits]),
caxis=xrtp.XYCAxis(
label=r"energy",
unit=r"eV"),
title=r"Source")
plots.append(plot00)

fwhmFormatStrF = '%.2f'
plotsFSM2 = []
for i, dq in enumerate(beamLine.fsm2.dqs):
    plot2 = xrtp.XYCPlot(
        'beamFSM2_{0:02d}'.format(i), (1,),
        xaxis=xrtp.XYCAxis(
            r'$x$', u'µm'),#, limits=[-xyLimits,xyLimits]
        yaxis=xrtp.XYCAxis(
            r'$z$', u'µm'),
        title=beamLine.fsm2.name+'_{0:02d}'.format(i),
        caxis=xrtp.XYCAxis('energy', 'eV',offset=E0,limits=[E0-2,E0+2])
        )
    plot2.xaxis.fwhmFormatStr = fwhmFormatStrF
    plot2.yaxis.fwhmFormatStr = fwhmFormatStrF
    plot2.textPanel = plot2.fig.text(
        0.2, 0.75, '', transform=plot2.fig.transFigure, size=14, color='r',
        ha='left')
    plot2.textPanelTemplate = '{0}: d$q=${1:+.0f} mm'.format('{0}', dq)
    plots.append(plot2)
    plotsFSM2.append(plot2)    
return plots,plotsFSM2

def plot_generator(plots, plotsFSM2, beamLine):

figDF = plt.figure(figsize=(7, 5), dpi=72)
ax1 = plt.subplot(111)
ax1.set_title(r'FWHM size of beam cross-section near focal position')
ax1.set_xlabel(r'd$q$ (mm)', fontsize=14)
ax1.set_ylabel(u'FWHM size (µm)', fontsize=14)
prefix = 'CRL-indiv-'
for plot in plots:
    fileName = '{0}{1}'.format(
        prefix, plot.title)
    plot.saveName = fileName + '.png'

plot.persistentName = fileName + '.pickle'

    try:
        plot.textPanel.set_text(
            plot.textPanelTemplate.format(material.elements[0].name))
    except AttributeError:
        pass
yield
for i, lens in enumerate(beamLine.lenses):
    True
           
qCurve = []
xCurve = []
zCurve = []
for dq, plot in zip(beamLine.fsm2.dqs, plotsFSM2):
    qCurve.append(dq)            
    xCurve.append(plot.dx)
    zCurve.append(plot.dy)
ax1.plot(qCurve, xCurve, 'o', label='Horizontal direction {0}, n={1:0d}'.format(
                material.elements[0].name, i+1))  
ax1.plot(qCurve, zCurve, '+', label='Vertical direction {0}, n={1:0d}'.format(
                material.elements[0].name, i+1))
ax1.set_xlabel(u'Shift (mm)', fontsize=14)
ax1.set_ylabel(u'FWHM size (mm)', fontsize=14)
plt.legend()          
plt.show()
plt.savefig('depthOfFocus2.png',dpi=200)

def main():
beamLine = build_beamline()
plots, plotsFSM2 = define_plots(beamLine)
xrtr.run_ray_tracing(
plots, repeats=1, updateEvery=1,generator=plot_generator,
generatorArgs=[plots, plotsFSM2, beamLine],
beamLine=beamLine)

if name == 'main':
main()

[yangfg-bsrf]

In this codes, the slope error is generated by the other fille and the data is loaded by the read_data(self) function, the error is added by modified the local_z1(self, x, y) and local_n1(self, x, y) functions.

Hello,
Where should I look at? Which part of the code is yours? Where in the code "the slope error is R=100um"? What is "CRL’s R" and where in the code you use it?

[kklmn]

I see you use persistentName parameter. Please make sure that you understand what it does and pay attention to the warning here. It might be the reason for the misbehavior of your code.

Please edit your posts so that the code parts are formatted properly.

[Faryzww]

Oh，there is a # before persistentName, I don't know why a lot of # are missing when I paste the codes here. The codes can run properly without a warning on my computer. But the result is wrong.

[Faryzww]

@goxrt @kklmn
Yes, but I don't know whether it is right or not, because the final reult is wrong, and I don't know why.

In this codes, the slope error is generated by the other fille and the data is loaded by the read_data(self) function, the error is added by modified the local_z1(self, x, y) and local_n1(self, x, y) functions.

[kklmn]

Please edit your post with tripple quotes, I sent you a link for explanations earlier.

Before you post your message please preview it first.

What are the reasons to believe that the files you load do really contain what you expect? Even if they do, how do you expect I could test it?

[Faryzww]

@kklmn
Sorry, maybe I didn't express well, and there's some misunderstanding between us. So I will change my question: how to add a slope error to ideal CRLs properly? Can you give me an example?

[kklmn]

Yes, you were right to assume that the mirror example should work for lenses.

What can be wrong in your implementation is that you add the same error on both sides of each lens. The error on the other side of a lens is flipped together with the normal. The two normals on both surfaces (local z directions) are anti-parallel, the local y directions are equal and the x's are again antiparallel. The resulting effect depends on your definition of the error tables in use.

If you want me to look at your scripts, please make them readable.

[Faryzww]

Yes, I add the same error on both sides of each lens. In my test code, the error is rotational symmetrical, I don't think the fipping will make the error added wrong. I will try to paste my script again. God bless me.
import os, sys;

sys.path.append(os.path.join('..', '..', '..'))
import numpy as np
from scipy import ndimage
import xrt.backends.raycing.oes as roe

def Read_distorted_surface(fname1, fname2, fname3):
x = np.loadtxt(fname1, unpack=True)
y = np.loadtxt(fname2, unpack=True)
z = np.loadtxt(fname3, unpack=True)
return x, y, z

class DoubleParaboloidLens_Distorted(roe.DoubleParaboloidLens):

hiddenMethods = roe.DoubleParaboloidLens.hiddenMethods + ['double_refract']  

def __init__(self, *args, **kwargs):
    kwargs = self.__pop_kwargs(**kwargs)
    roe.DoubleParaboloidLens.__init__(self, *args, **kwargs)
    self.read_data()

def __pop_kwargs(self, **kwargs):
    self.fname1 = kwargs.pop('fname1') 
    self.fname2 = kwargs.pop('fname2')
    self.fname3 = kwargs.pop('fname3')
    self.get_distorted_surface = kwargs.pop('get_distorted_surface')
    return kwargs

def assign_auto_material_kind(self, material):
    material.kind = 'lens'

def read_data(self):
    if self.get_distorted_surface is None:
        return
    self.warpX, self.warpY, self.warpZ = Read_distorted_surface(self.fname1, self.fname2, self.fname3)
    self.warpNX, self.warpNY = len(self.warpX), len(self.warpY)
    self.limPhysX0 = np.min(self.warpX), np.max(self.warpX)
    self.limPhysY0 = np.min(self.warpY), np.max(self.warpY)
    self.get_surface_limits()
    self.warpA, self.warpB = np.gradient(self.warpZ)
    dx = self.warpX[1] - self.warpX[0]
    dy = self.warpY[1] - self.warpY[0]
    self.warpA = np.arctan(self.warpA / dx)
    self.warpB = np.arctan(self.warpB / dy)

    self.warpSplineZ = ndimage.spline_filter(self.warpZ)
    self.warpSplineA = ndimage.spline_filter(self.warpA)
    self.warpSplineB = ndimage.spline_filter(self.warpB)

def local_z1(self, x, y):

    z = (x ** 2 + y ** 2) / (4 * self.focus)
    coords = np.array(
        [(x - self.limPhysX0[0]) /
         (self.limPhysX0[1] - self.limPhysX0[0]) * (self.warpNX - 1),
         (y - self.limPhysY[0]) /
         (self.limPhysY0[1] - self.limPhysY0[0]) * (self.warpNY - 1)])

    z = z + ndimage.map_coordinates(self.warpSplineZ, coords, prefilter=True)
    if self.zmax is not None:
        z[z > self.zmax] = self.zmax
    return z

def local_z2(self, x, y):
    return self.local_z1(x, y)


def local_n1(self, x, y):
    Ax = -self.warpSplineA
    By = -self.warpSplineB
    coords = np.array(
        [(x - self.limPhysX0[0]) /
         (self.limPhysX0[1] - self.limPhysX0[0]) * (self.warpNX - 1),
         (y - self.limPhysY0[0]) /
         (self.limPhysY0[1] - self.limPhysY0[0]) * (self.warpNY - 1)])
    a0 = ndimage.map_coordinates(Ax, coords, prefilter=True)
    b0 = ndimage.map_coordinates(By, coords, prefilter=True)
    if self.zmax is not None:
        z = (x ** 2 + y ** 2) / (4 * self.focus)
        z = z + ndimage.map_coordinates(self.warpSplineZ, coords, prefilter=True)

        a = -x / (2 * self.focus) + a0 
        b = -y / (2 * self.focus) + b0 

        if isinstance(a, np.ndarray):
            a[z > self.zmax] = 0
        if isinstance(b, np.ndarray):
            b[z > self.zmax] = 0

        c = np.ones_like(x)
        norm = (a ** 2 + b ** 2 + 1) ** 0.5
    return [a / norm, b / norm, c / norm]

def local_n2(self, x, y):
    return self.local_n1(x, y)

[Faryzww]

import os, sys; sys.path.append(os.path.join('..', '..', '..'))
import os, sys; sys.path.append(os.path.join('C:\Users\Lenovo\Desktop\XRT_Simulation\surface error\Error_Generation'))
import numpy as np
import matplotlib.pyplot as plt
import xrt.backends.raycing as raycing
import xrt.backends.raycing.sources as rs
import xrt.backends.raycing.oes as roe
import xrt.backends.raycing.apertures as ra
import xrt.backends.raycing.run as rr
import xrt.backends.raycing.materials as rm
import xrt.backends.raycing.screens as rsc
import xrt.plotter as xrtp
import xrt.runner as xrtr
import Distorted_Element as Distorted

Lens = roe.DoubleParaboloidLens
Lens_Err=Distorted.DoubleParaboloidLens_Distorted

R = 0.1
focus = R/2
ztatal = 0.255
t=0.03
zmax = (ztatal - t)/2
dz = 3.
E0 = 12400.
f = 1770.
p = 50000
print('q = {0:.2f}'.format(1/(1/f-1/p)))
xyLimits = -5, 5

mBeryllium = rm.Material('Be', rho=1.845, kind='lens')
material = mBeryllium

dE = E0*0.0001/2
eMin0,eMax0 = E0-dE,E0+dE
nrays = 1e4
kwargs_B4 = dict(
nrays=nrays, center=[0, 0, 0],
eE=6.0, eI=0.2, eEspread=0.00111,
eEpsilonX=0.02728, eEpsilonZ=0.00276,
betaX=10.12, betaZ=9.64,
xPrimeMax=0.0012, zPrimeMax=0.0012,
targetE=[E0, 3], eMin=eMin0, eMax=eMax0,
uniformRayDensity=True,
K=r"auto",period=32,n=156)

def build_beamline(nrays=1e4):

beamLine = raycing.BeamLine()    
beamLine.source = rs.Undulator(beamLine,**kwargs_B4) 
beamLine.fsm1 = rsc.Screen(beamLine, 'FSM1', (0, p - 500, 0))        
beamLine.roundAperture01 = ra.RoundAperture(
    bl=beamLine,
    name=None,
    center=[0, 49999, 0],
    r=0.06)
beamLine.lenses = []     
lens_error = 1    
Surface_name='CRL_try_'   
kwargs_OE1 = dict(pitch=np.pi/2, t=t, material=material,
                limPhysX=[-0.12/2, 0.12/2], limPhysY=[-0.12/2, 0.12/2], shape='round',
                focus=focus, zmax=zmax, alarmLevel=0.1)      
ilens = 0
while True: 
    center=[0, p + dz*ilens, 0]           
    if lens_error==0:
        lens = Lens(bl=beamLine,name = 'Lens{0:02d}'.format(ilens), 
                center=center,**kwargs_OE1)           
    else:
        kwargs_OE1['get_distorted_surface'] ='error'
        kwargs_OE1['fname1'] =Surface_name+'X.txt'
        kwargs_OE1['fname2'] =Surface_name+'Y.txt'
        kwargs_OE1['fname3'] =Surface_name+'Z.txt'              
        lens = Lens_Err(bl=beamLine,name ='Lens{0:02d}'.format(ilens), 
                center=center,**kwargs_OE1)                  
    beamLine.lenses.append(lens)
    if ilens == 0:
        nCRL = lens.get_nCRL(f, E0)           
    ilens += 1
    if nCRL - ilens < 0.5:  
        break 

beamLine.fsm2 = rsc.Screen(beamLine,'FSM2')
beamLine.fsm2.dqs = np.linspace(-1800, 100, 11)

return beamLine

def run_process(beamLine):

beamSource = beamLine.source.shine()
outDict = {'beamSource': beamSource} 

beamFSM1 = beamLine.fsm1.expose(beamSource)
outDict['beamFSM1'] = beamFSM1

beamIn = beamSource
for ilens, lens in enumerate(beamLine.lenses):
    lglobal, llocal1, llocal2 = lens.double_refract(beamIn, needLocal=True)
    beamIn = lglobal
    strl = '_{0:02d}'.format(ilens)
    outDict['beamLensGlobal'+strl] = lglobal
    outDict['beamLensLocal1'+strl] = llocal1
    outDict['beamLensLocal2'+strl] = llocal2
    
if Lens==roe.DoubleParaboloidLens:
    nFactor = 0.5
elif Lens==roe.ParabolicCylinderFlatLens:
    nFactor = 2.
else:
    nFactor = 1.
delta=(1. - material.get_refractive_index(E0).real)     
f_real= 2 * focus / len(beamLine.lenses) /delta*nFactor    
q=1/(1/f_real-1/p)
print('Focal length of the lenses ={0:3} mm'.format(f_real))
print('The number of the lenses ={0:3} '.format(ilens))
print('The number of the lenses ={0:3} '.format(ilens))


for i, dq in enumerate(beamLine.fsm2.dqs):        
    beamLine.fsm2.center[1] = p + q + dq
    outDict['beamFSM2_{0:02d}'.format(i)] = beamLine.fsm2.expose(lglobal)
          
return outDict

rr.run_process = run_process

def define_plots(beamLine):
plots = []
xrtp.yTextPosNraysR = 0.82
xrtp.yTextPosNrays1 = 0.52
xyLimits = 50
plot00 = xrtp.XYCPlot(
beam=r"beamFSM1",
xaxis=xrtp.XYCAxis(
label=r"x",unit=r"$\mu$m",limits=[-xyLimits,xyLimits]),
yaxis=xrtp.XYCAxis(
label=r"z",unit=r"$\mu$m",limits=[-xyLimits,xyLimits]),
caxis=xrtp.XYCAxis(
label=r"energy",
unit=r"eV"),
title=r"Source")
plots.append(plot00)

fwhmFormatStrF = '%.2f'
plotsFSM2 = []
for i, dq in enumerate(beamLine.fsm2.dqs):
    plot2 = xrtp.XYCPlot(
        'beamFSM2_{0:02d}'.format(i), (1,),
        xaxis=xrtp.XYCAxis(
            r'$x$', u'µm'),
        yaxis=xrtp.XYCAxis(
            r'$z$', u'µm'),
        title=beamLine.fsm2.name+'_{0:02d}'.format(i),
        caxis=xrtp.XYCAxis('energy', 'eV',offset=E0,limits=[E0-2,E0+2])
        )
    plot2.xaxis.fwhmFormatStr = fwhmFormatStrF
    plot2.yaxis.fwhmFormatStr = fwhmFormatStrF
    plot2.textPanel = plot2.fig.text(
        0.2, 0.75, '', transform=plot2.fig.transFigure, size=14, color='r',
        ha='left')
    plot2.textPanelTemplate = '{0}: d$q=${1:+.0f} mm'.format('{0}', dq)
    plots.append(plot2)
    plotsFSM2.append(plot2)    
return plots,plotsFSM2

def plot_generator(plots, plotsFSM2, beamLine):

ax1 = plt.subplot(111)
ax1.set_title(r'FWHM size of beam cross-section near focal position')
ax1.set_xlabel(r'd$q$ (mm)', fontsize=14)
ax1.set_ylabel(u'FWHM size (µm)', fontsize=14)
prefix = 'CRL-indiv-'
for plot in plots:
    fileName = '{0}{1}'.format(
        prefix, plot.title)
    plot.saveName = fileName + '.png'
    try:
        plot.textPanel.set_text(
            plot.textPanelTemplate.format(material.elements[0].name))
    except AttributeError:
        pass
yield
for i, lens in enumerate(beamLine.lenses):
    True
           
qCurve = []
xCurve = []
zCurve = []
for dq, plot in zip(beamLine.fsm2.dqs, plotsFSM2):
    qCurve.append(dq)            
    xCurve.append(plot.dx)
    zCurve.append(plot.dy)
ax1.plot(qCurve, xCurve, 'o', label='Horizontal direction {0}, n={1:0d}'.format(
                material.elements[0].name, i+1))  
ax1.plot(qCurve, zCurve, '+', label='Vertical direction {0}, n={1:0d}'.format(
                material.elements[0].name, i+1))
ax1.set_xlabel(u'Shift (mm)', fontsize=14)
ax1.set_ylabel(u'FWHM size (mm)', fontsize=14)
plt.legend()          
plt.show()
plt.savefig('depthOfFocus2.png',dpi=200)

def main():
beamLine = build_beamline()
plots, plotsFSM2 = define_plots(beamLine)
xrtr.run_ray_tracing(
plots, repeats=1, updateEvery=1,generator=plot_generator,
generatorArgs=[plots, plotsFSM2, beamLine],
beamLine=beamLine)

if name == 'main':
main()