search

kklmn / xrt

Package xrt (XRayTracer) is a python software library for ray tracing and wave propagation in x-ray regime. It is primarily meant for modeling synchrotron sources, beamlines and beamline elements.
MIT License
83 stars 30 forks source link

Adding slope errors to CRLs (double paraboloid lens) #40

Open Faryzww opened 5 years ago

Faryzww commented 5 years ago

Hi, I have tried to add slope errors to CRLs, and assumed every surface had the same slope error. The slope error is added according to mirrors (see the following codes). In order to know whether it works right, I add the slope error to change the CRL’s R to another value (for example, R of CRLs is 100um, the slope error is R=100um, and then the final R will be 50um). But unfortunately, the result is absolutely wrong.

------------------------------------------------------------

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')  # 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):
    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)
Faryzww commented 5 years ago

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 commented 5 years ago

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 commented 5 years ago

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 commented 5 years ago

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 commented 5 years ago

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 commented 5 years ago

In this codes, the slope error is generated by the other fille and the data is loaded by the _readdata(self) function, the error is added by modified the _localz1(self, x, y) and _localn1(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 commented 5 years ago

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 commented 5 years ago

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 commented 5 years ago

@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 commented 5 years ago

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 commented 5 years ago

@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 commented 5 years ago

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 commented 5 years ago

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 commented 5 years ago

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()