From 02a6215552326723e23cbca735b13101d02e69ca Mon Sep 17 00:00:00 2001 From: Roman Chernikov Date: Mon, 22 Jun 2026 21:55:25 -0400 Subject: [PATCH 1/6] Basic mirror reflectivity calculation with xraydb --- env/python/requirements-optional.txt | 3 +- .../examples/SRWLIB_Example09-materials.py | 346 +++++++++++++++ env/python/srwpy/uti_mtrl.py | 413 ++++++++++++++++++ 3 files changed, 761 insertions(+), 1 deletion(-) create mode 100644 env/python/srwpy/examples/SRWLIB_Example09-materials.py create mode 100644 env/python/srwpy/uti_mtrl.py diff --git a/env/python/requirements-optional.txt b/env/python/requirements-optional.txt index 3e01382d..ea211af2 100644 --- a/env/python/requirements-optional.txt +++ b/env/python/requirements-optional.txt @@ -3,4 +3,5 @@ oasys_srw primme pykern mpld3 -scikit-image \ No newline at end of file +scikit-image +xraydb diff --git a/env/python/srwpy/examples/SRWLIB_Example09-materials.py b/env/python/srwpy/examples/SRWLIB_Example09-materials.py new file mode 100644 index 00000000..894afccb --- /dev/null +++ b/env/python/srwpy/examples/SRWLIB_Example09-materials.py @@ -0,0 +1,346 @@ +############################################################################# +# SRWLIB Example#9_Materials: Simulating propagation of a Gaussian X-ray beam +# through a Beamline Containing Mirrors made of different materials +# v 0.04 +############################################################################# + +from __future__ import print_function #Python 2.7 compatibility + +try: #OC15112022 + import sys + sys.path.append('../') + from srwlib import * + from uti_plot import * + from uti_mtrl import * +except: + from srwpy.srwlib import * + from srwpy.uti_plot import * + from srwpy.uti_mtrl import * +import os +import copy + +print('SRWLIB Python Example # 9 with material properties:') +print('Simulating propagation of a coherent Gaussian X-ray beam through Pt and B4C mirrors') + + + +#******************* Reflectivity Parameters + +_n_ph_en=100 +_n_ang=1000 +_n_comp=2 +_ph_en_start=12000 +_ph_en_fin=13000 +_ph_en_scale_type="lin" +_ang_start=0.001 +_ang_fin=0.01 +_ang_scale_type="lin" +avgEn = 12400 + + +Pt_refl=calc_refl_arr('Pt', _n_ang=_n_ang, _n_ph_en=_n_ph_en, _n_comp=_n_comp, _ph_en_start=_ph_en_start, _ph_en_fin=_ph_en_fin, _ph_en_scale_type=_ph_en_scale_type, + _ang_start=_ang_start, _ang_fin=_ang_fin, _ang_scale_type=_ang_scale_type) +B4C_refl=calc_refl_arr('B4C', _dens=2.5, _n_ang=_n_ang, _n_ph_en=_n_ph_en, _n_comp=_n_comp, _ph_en_start=_ph_en_start, _ph_en_fin=_ph_en_fin, _ph_en_scale_type=_ph_en_scale_type, + _ang_start=_ang_start, _ang_fin=_ang_fin, _ang_scale_type=_ang_scale_type) + + + +#**********************Input Parameters and Structures +#***********Folder and Data File Names +strDataFolderName = 'data_example_09-materials' #data sub-folder name +strIntOutFileName = 'ex09_res_int_in.dat' #initial wavefront intensity distribution output file name +strPhOutFileName = 'ex09_res_phase_in.dat' #initial wavefront phase output file name +strIntPropOutFileName = 'ex09_res_int_prop.dat' #propagated wavefront intensity distribution output file name +strPhasePropOutFileName = 'ex09_res_phase_prop.dat' #propagated wavefront phase output file name +strIntPropOutFileName2 = 'ex09_res_int_prop_Pt.dat' #propagated wavefront intensity distribution output file name +strPhPropOutFileName2 = 'ex09_res_phase_prop_Pt.dat' #propagated wavefront phase output file name +strIntPropOutFileName3 = 'ex09_res_int_prop_B4C.dat' #propagated wavefront intensity distribution output file name +strPhPropOutFileName3 = 'ex09_res_phase_prop_B4C.dat' #propagated wavefront phase output file name + + +#***********Gaussian Beam Source +GsnBm = SRWLGsnBm() #Gaussian Beam structure (just parameters) +GsnBm.x = 0 #Transverse Positions of Gaussian Beam Center at Waist [m] +GsnBm.y = 0 +GsnBm.z = 0 #Longitudinal Position of Waist [m] +GsnBm.xp = 0 #Average Angles of Gaussian Beam at Waist [rad] +GsnBm.yp = 0 +GsnBm.avgPhotEn = avgEn #Photon Energy [eV] +GsnBm.pulseEn = 0.001 #Energy per Pulse [J] - to be corrected +GsnBm.repRate = 1 #Rep. Rate [Hz] - to be corrected +GsnBm.polar = 1 #1- linear horizontal +GsnBm.sigX = 23e-06/2.35 #Horiz. RMS size at Waist [m] +GsnBm.sigY = GsnBm.sigX #Vert. RMS size at Waist [m] + +constConvRad = 1.23984186e-06/(4*3.1415926536) +rmsAngDiv = constConvRad/(GsnBm.avgPhotEn*GsnBm.sigX) #RMS angular divergence [rad] +print('RMS Source Size:', round(GsnBm.sigX*1.e+06, 3), 'micro-m; RMS Divergence:', round(rmsAngDiv*1.e+06, 3), 'micro-rad') + +GsnBm.sigT = 10e-15 #Pulse duration [fs] (not used?) +GsnBm.mx = 0 #Transverse Gauss-Hermite Mode Orders +GsnBm.my = 0 + +#***********Initial Wavefront +wfr = SRWLWfr() #Initial Electric Field Wavefront +wfr.allocate(1, 100, 100) #Numbers of points vs Photon Energy (1), Horizontal and Vertical Positions (dummy) +wfr.mesh.zStart = 25 #Longitudinal Position [m] at which initial Electric Field has to be calculated, i.e. the position of the first optical element +wfr.mesh.eStart = GsnBm.avgPhotEn #Initial Photon Energy [eV] +wfr.mesh.eFin = GsnBm.avgPhotEn #Final Photon Energy [eV] + +wfr.unitElFld = 2 #Electric field units: 0- arbitrary, 1- sqrt(Phot/s/0.1%bw/mm^2), 2- sqrt(J/eV/mm^2) or sqrt(W/mm^2), depending on representation (freq. or time) + +distSrc_VFM = wfr.mesh.zStart - GsnBm.z +#Horizontal and Vertical Position Range for the Initial Wavefront calculation +#can be used to simulate the First Aperture +firstHorAp = 8.*rmsAngDiv*distSrc_VFM #[m] +firstVertAp = firstHorAp #[m] + +wfr.mesh.xStart = -0.5*firstHorAp #Initial Horizontal Position [m] +wfr.mesh.xFin = 0.5*firstHorAp #Final Horizontal Position [m] +wfr.mesh.yStart = -0.5*firstVertAp #Initial Vertical Position [m] +wfr.mesh.yFin = 0.5*firstVertAp #Final Vertical Position [m]s + +sampFactNxNyForProp = 4 #sampling factor for adjusting nx, ny (effective if > 0) +arPrecPar = [sampFactNxNyForProp] + +wfr.partBeam.partStatMom1.x = GsnBm.x #Some information about the source in the Wavefront structure +wfr.partBeam.partStatMom1.y = GsnBm.y +wfr.partBeam.partStatMom1.z = GsnBm.z +wfr.partBeam.partStatMom1.xp = GsnBm.xp +wfr.partBeam.partStatMom1.yp = GsnBm.yp + + +#**********************Calculating Initial Wavefront and extracting Intensity: +srwl.CalcElecFieldGaussian(wfr, GsnBm, arPrecPar) +arI0 = array('f', [0]*wfr.mesh.nx*wfr.mesh.ny) #"flat" array to take 2D intensity data +srwl.CalcIntFromElecField(arI0, wfr, 6, 0, 3, wfr.mesh.eStart, 0, 0) #extracts intensity +srwl_uti_save_intens_ascii(arI0, wfr.mesh, os.path.join(os.getcwd(), strDataFolderName, strIntOutFileName), 0) + +arI0x = array('f', [0]*wfr.mesh.nx) #array to take 1D intensity data +srwl.CalcIntFromElecField(arI0x, wfr, 6, 0, 1, wfr.mesh.eStart, 0, 0) #extracts intensity + +arP0 = array('d', [0]*wfr.mesh.nx*wfr.mesh.ny) #"flat" array to take 2D phase data (note it should be 'd') +srwl.CalcIntFromElecField(arP0, wfr, 0, 4, 3, wfr.mesh.eStart, 0, 0) #extracts radiation phase +srwl_uti_save_intens_ascii(arP0, wfr.mesh, os.path.join(os.getcwd(), strDataFolderName, strPhOutFileName), 0, ['', 'Horizontal Position', 'Vertical Position', 'Phase'], _arUnits=['', 'm', 'm', 'rad']) + +mesh0 = deepcopy(wfr.mesh) + + + +#***********Optical Elements and Propagation Parameters +#Sequence of Optical Elements: +# +# +# +# +# +# + + + +distVFM_HFM = 1 #Distance from VFM to HFM [m] +distHFM_Samp = 1 #Distance from HFM to Sample [m] + + +lenKB = 0.1 #Length of VFM and HFM (same for each) [m] +angKB = 0.004 #b4c_ca #2.521e-3 #Incident Angle of VFM and HFM [rad] + + + +#Aperture of KB system +opApKB = SRWLOptA('r', 'a', lenKB*angKB, lenKB*angKB) + + +#VFM simulated by Ideal Lens: +#opVFM = SRWLOptL(_Fy=(distSrc_M1 + distM1_VFM)*(distVFM_HFM + distHFM_Samp)/(distSrc_M1 + distM1_VFM + distVFM_HFM + distHFM_Samp)) +#VFM simulated by Extended Elliptical Mirror: +opVFM = SRWLOptMirEl(_p=(distSrc_VFM), _q=(distVFM_HFM + distHFM_Samp), _ang_graz=angKB, _size_tang=lenKB, _size_sag=10.e-03, + _nvx=0, _nvy=cos(angKB), _nvz=-sin(angKB), _tvx=0, _tvy=-sin(angKB)) + +#Drift from VFM to HFM +opDrVFM_HFM = SRWLOptD(distVFM_HFM) + +#VFM simulated by Extended Elliptical Mirror: +opHFM = SRWLOptMirEl(_p=(distSrc_VFM + distVFM_HFM), _q=distHFM_Samp, _ang_graz=angKB, _size_tang=lenKB, _size_sag=10.e-03, + _nvx=cos(angKB), _nvy=0, _nvz=-sin(angKB), _tvx=-sin(angKB), _tvy=0) + +#Drift from HFM to Sample +opDrHFM_Samp = SRWLOptD(distHFM_Samp) + + +#Wavefront Propagation Parameters: +#[0]: Auto-Resize (1) or not (0) Before propagation +#[1]: Auto-Resize (1) or not (0) After propagation +#[2]: Relative Precision for propagation with Auto-Resizing (1. is nominal) +#[3]: Allow (1) or not (0) for semi-analytical treatment of the quadratic (leading) phase terms at the propagation +#[4]: Do any Resizing on Fourier side, using FFT, (1) or not (0) +#[5]: Horizontal Range modification factor at Resizing (1. means no modification) +#[6]: Horizontal Resolution modification factor at Resizing +#[7]: Vertical Range modification factor at Resizing +#[8]: Vertical Resolution modification factor at Resizing +#[9]: Type of wavefront Shift before Resizing (not yet implemented) +#[10]: New Horizontal wavefront Center position after Shift (not yet implemented) +#[11]: New Vertical wavefront Center position after Shift (not yet implemented) +# [0][1][2] [3][4] [5] [6] [7] [8] [9][10][11] +#prParInit = [0, 0, 1., 1, 0, 2., 5., 2., 3., 0, 0, 0] +prParInit = [0, 0, 1., 1, 0, 2., 2., 2., 2., 0, 0, 0] +prPar0 = [0, 0, 1., 1, 0, 1., 2, 1., 2, 0, 0, 0] +prPar1 = [0, 0, 1., 1, 0, 1, 2, 0.1, 2, 0, 0, 0] +#prParPost = [0, 0, 1., 1, 0, 0.06, 3., 0.1, 2., 0, 0, 0] + +#NOTE: in this case of simulation, it can be enough to define the precision parameters only Before and After +#the propagation through the entire Beamline. However, if necessary, different propagation parameters can be specified +#for each optical element. +#The optimal values of propagation parameters may depend on photon energy and optical layout. + +#"Beamline" - a sequenced Container of Optical Elements (together with the corresponding wavefront propagation parameters, +#and the "post-propagation" wavefront resizing parameters for better viewing). +# optBL = SRWLOptC([opVFM, opTrErVFM, opDrVFM_HFM, opHFM, opTrErHFM, opDrHFM_Samp], +# [prPar0, prPar0, prPar0, prPar0, prPar0, prParPost]) +optBL = SRWLOptC([opApKB, opVFM, opDrVFM_HFM, opHFM, opDrHFM_Samp], + [prParInit, prPar0, prPar0, prPar0, prPar1]) + + +# # Beamline with material properties +# #VFM simulated by Ideal Lens: +# #opVFM = SRWLOptL(_Fy=(distSrc_M1 + distM1_VFM)*(distVFM_HFM + distHFM_Samp)/(distSrc_M1 + distM1_VFM + distVFM_HFM + distHFM_Samp)) +# #VFM simulated by Extended Elliptical Mirror: +opVFM_Pt = SRWLOptMirEl(_p=(distSrc_VFM), _q=(distVFM_HFM + distHFM_Samp), _ang_graz=angKB, _size_tang=lenKB, _size_sag=10.e-03, + _nvx=0, _nvy=cos(angKB), _nvz=-sin(angKB), _tvx=0, _tvy=-sin(angKB), + _refl=Pt_refl, _n_ang=_n_ang, _n_ph_en=_n_ph_en, _n_comp=_n_comp, _ph_en_start=_ph_en_start, _ph_en_fin=_ph_en_fin, _ph_en_scale_type=_ph_en_scale_type, + _ang_start=_ang_start, _ang_fin=_ang_fin, _ang_scale_type=_ang_scale_type) +# # opTrErVFM_Pt = srwl_opt_setup_surf_height_1d(heightProfData, _dim='y', _ang=angKB, _amp_coef=1) + +# #HFM simulated by Ideal Lens: +# #opHFM = SRWLOptL(_Fx=(distSrc_M1 + distM1_VFM + distVFM_HFM)*distHFM_Samp/(distSrc_M1 + distM1_VFM + distVFM_HFM + distHFM_Samp)) +# #VFM simulated by Extended Elliptical Mirror: +opHFM_Pt = SRWLOptMirEl(_p=(distSrc_VFM + distVFM_HFM), _q=distHFM_Samp, _ang_graz=angKB, _size_tang=lenKB, _size_sag=10.e-03, + _nvx=cos(angKB), _nvy=0, _nvz=-sin(angKB), _tvx=-sin(angKB), _tvy=0, + _refl=Pt_refl, _n_ang=_n_ang, _n_ph_en=_n_ph_en, _n_comp=_n_comp, _ph_en_start=_ph_en_start, _ph_en_fin=_ph_en_fin, _ph_en_scale_type=_ph_en_scale_type, + _ang_start=_ang_start, _ang_fin=_ang_fin, _ang_scale_type=_ang_scale_type) +# opTrErHFM_Pt = srwl_opt_setup_surf_height_1d(heightProfData, _dim='x', _ang=angKB, _amp_coef=1) + +# #"Beamline" - a sequenced Container of Optical Elements (together with the corresponding wavefront propagation parameters, +# #and the "post-propagation" wavefront resizing parameters for better viewing). +# # optBL_Pt = SRWLOptC([opVFM_Pt, opTrErVFM_Pt, opDrVFM_HFM, opHFM_Pt, opTrErHFM_Pt, opDrHFM_Samp], +# # [prPar0, prPar0, prPar0, prPar0, prPar0, prParPost]) +optBL_Pt = SRWLOptC([opApKB, opVFM_Pt, opDrVFM_HFM, opHFM_Pt, opDrHFM_Samp], + [prParInit, prPar0, prPar0, prPar0, prPar1]) + + + +# Boron Carbide Mirrors +opVFM_B4C = SRWLOptMirEl(_p=(distSrc_VFM), _q=(distVFM_HFM + distHFM_Samp), _ang_graz=angKB, _size_tang=lenKB, _size_sag=10.e-03, + _nvx=0, _nvy=cos(angKB), _nvz=-sin(angKB), _tvx=0, _tvy=-sin(angKB), + _refl=B4C_refl, _n_ang=_n_ang, _n_ph_en=_n_ph_en, _n_comp=_n_comp, _ph_en_start=_ph_en_start, _ph_en_fin=_ph_en_fin, _ph_en_scale_type=_ph_en_scale_type, + _ang_start=_ang_start, _ang_fin=_ang_fin, _ang_scale_type=_ang_scale_type) + +opHFM_B4C = SRWLOptMirEl(_p=(distSrc_VFM + distVFM_HFM), _q=distHFM_Samp, _ang_graz=angKB, _size_tang=lenKB, _size_sag=10.e-03, + _nvx=cos(angKB), _nvy=0, _nvz=-sin(angKB), _tvx=-sin(angKB), _tvy=0, + _refl=B4C_refl, _n_ang=_n_ang, _n_ph_en=_n_ph_en, _n_comp=_n_comp, _ph_en_start=_ph_en_start, _ph_en_fin=_ph_en_fin, _ph_en_scale_type=_ph_en_scale_type, + _ang_start=_ang_start, _ang_fin=_ang_fin, _ang_scale_type=_ang_scale_type) + +optBL_B4C = SRWLOptC([opApKB, opVFM_B4C, opDrVFM_HFM, opHFM_B4C, opDrHFM_Samp], + [prParInit, prPar0, prPar0, prPar0, prPar1]) + + + + +#***********Wavefront Propagation +arI1 = None; arI1x = None; arI1y = None +arP1 = None +mesh1 = None + + +#Duplicating wavefront (by re-creating all objects/arrays): +wfr1 = copy.deepcopy(wfr) +wfr2 = copy.deepcopy(wfr) +wfr3 = copy.deepcopy(wfr) + + + +print(' Propagating Wavefront (through two elliptical mirrors)... ', end='') +srwl.PropagElecField(wfr1, optBL) +print('done') +print(' Saving resulting Intensity and Phase data to files ... ', end='') +mesh1 = deepcopy(wfr1.mesh) +arI1 = array('f', [0]*mesh1.nx*mesh1.ny) #"flat" array to take 2D intensity data +srwl.CalcIntFromElecField(arI1, wfr1, 6, 0, 3, mesh1.eStart, 0, 0) #extracts intensity +srwl_uti_save_intens_ascii(arI1, mesh1, os.path.join(os.getcwd(), strDataFolderName, strIntPropOutFileName), 0) + +arI1x = array('f', [0]*mesh1.nx) #array to take 1D intensity data +srwl.CalcIntFromElecField(arI1x, wfr1, 6, 0, 1, mesh1.eStart, 0, 0) #extracts intensity +arI1y = array('f', [0]*mesh1.ny) #array to take 1D intensity data +srwl.CalcIntFromElecField(arI1y, wfr1, 6, 0, 2, mesh1.eStart, 0, 0) #extracts intensity +arP1 = array('d', [0]*mesh1.nx*mesh1.ny) #"flat" array to take 2D phase data (note it should be 'd') +srwl.CalcIntFromElecField(arP1, wfr1, 0, 4, 3, mesh1.eStart, 0, 0) #extracts radiation phase +srwl_uti_save_intens_ascii(arP1, mesh1, os.path.join(os.getcwd(), strDataFolderName, strPhasePropOutFileName), 0, ['', 'Horizontal Position', 'Vertical Position', 'Phase'], _arUnits=['', 'm', 'm', 'rad']) +del wfr1 +print('done') + + +print(' Propagating Wavefront (through 2 Platinum Mirrors) ... ', end='') +srwl.PropagElecField(wfr2, optBL_Pt) +print('done') +print(' Saving resulting Intensity and Phase data to files ... ', end='') +mesh2 = deepcopy(wfr2.mesh) +arI2 = array('f', [0]*mesh2.nx*mesh2.ny) #"flat" array to take 2D intensity data +srwl.CalcIntFromElecField(arI2, wfr2, 6, 0, 3, mesh2.eStart, 0, 0) #extracts intensity +srwl_uti_save_intens_ascii(arI2, mesh2, os.path.join(os.getcwd(), strDataFolderName, strIntPropOutFileName2), 0) + +arI2x = array('f', [0]*mesh2.nx) #array to take 1D intensity data +srwl.CalcIntFromElecField(arI2x, wfr2, 6, 0, 1, mesh2.eStart, 0, 0) #extracts intensity +arI2y = array('f', [0]*mesh2.ny) #array to take 1D intensity data +srwl.CalcIntFromElecField(arI2y, wfr2, 6, 0, 2, mesh2.eStart, 0, 0) #extracts intensity +arP2 = array('d', [0]*mesh2.nx*mesh2.ny) #"flat" array to take 2D phase data (note it should be 'd') +srwl.CalcIntFromElecField(arP2, wfr2, 0, 4, 3, mesh2.eStart, 0, 0) #extracts radiation phase +srwl_uti_save_intens_ascii(arP2, mesh2, os.path.join(os.getcwd(), strDataFolderName, strPhPropOutFileName2), 0, ['', 'Horizontal Position', 'Vertical Position', 'Phase'], _arUnits=['', 'm', 'm', 'rad']) +del wfr2 +print('done') + +print(' Propagating Wavefront (through 2 B4C Mirrors) ... ', end='') +srwl.PropagElecField(wfr3, optBL_B4C) +print('done') +print(' Saving resulting Intensity and Phase data to files ... ', end='') +mesh3 = deepcopy(wfr3.mesh) +arI3 = array('f', [0]*mesh3.nx*mesh3.ny) +srwl.CalcIntFromElecField(arI3, wfr3, 6, 0, 3, mesh3.eStart, 0, 0) +srwl_uti_save_intens_ascii(arI3, mesh3, os.path.join(os.getcwd(), strDataFolderName, strIntPropOutFileName3), 0) + +arP3 = array('d', [0]*mesh3.nx*mesh3.ny) +srwl.CalcIntFromElecField(arP3, wfr3, 0, 4, 3, mesh3.eStart, 0, 0) +srwl_uti_save_intens_ascii(arP3, mesh3, os.path.join(os.getcwd(), strDataFolderName, strPhPropOutFileName3), 0, ['', 'Horizontal Position', 'Vertical Position', 'Phase'], _arUnits=['', 'm', 'm', 'rad']) +del wfr3 +print('done') + +#**********************Plotting results (requires 3rd party graphics package) +print(' Plotting the results (blocks script execution; close any graph windows to proceed) ... ', end='') +plotMesh0x = [1000*mesh0.xStart, 1000*mesh0.xFin, mesh0.nx] +plotMesh0y = [1000*mesh0.yStart, 1000*mesh0.yFin, mesh0.ny] +uti_plot2d1d(arI0, plotMesh0x, plotMesh0y, labels=['Horizontal Position [mm]', 'Vertical Position [mm]', 'Intensity Before Propagation']) + + +plotMesh1x = [1e+06*mesh1.xStart, 1e+06*mesh1.xFin, mesh1.nx] +plotMesh1y = [1e+06*mesh1.yStart, 1e+06*mesh1.yFin, mesh1.ny] +uti_plot2d1d(arI1, plotMesh1x, plotMesh1y, labels=['Horizontal Position [microns]', 'Vertical Position [microns]', 'Intensity After Propagation through 2 Elliptical Mirrors']) + +plotMesh2x = [1e+06*mesh2.xStart, 1e+06*mesh2.xFin, mesh2.nx] +plotMesh2y = [1e+06*mesh2.yStart, 1e+06*mesh2.yFin, mesh2.ny] +uti_plot2d1d(arI2, plotMesh2x, plotMesh2y, labels=['Horizontal Position [microns]', 'Vertical Position [microns]', 'Intensity After Propagation through 2 Platinum Elliptical Mirrors']) + +plotMesh3x = [1e+06*mesh3.xStart, 1e+06*mesh3.xFin, mesh3.nx] +plotMesh3y = [1e+06*mesh3.yStart, 1e+06*mesh3.yFin, mesh3.ny] +uti_plot2d1d(arI3, plotMesh3x, plotMesh3y, labels=['Horizontal Position [microns]', 'Vertical Position [microns]', 'Intensity After Propagation through 2 Boron Carbide Elliptical Mirrors']) + + + +# Calculate relative intensity after propagation through the material mirrors +print() +print("maximum intensity for perfect reflectivity: ", max(arI1)) +print("maximum intensity for Pt reflectivity: ", max(arI2)) +print("maximum intensity for B4C reflectivity: ", max(arI3)) +print("simulated relative intensity for Pt mirrors: ", max(arI2)/max(arI1)) +print("simulated relative intensity for B4C mirrors: ", max(arI3)/max(arI1)) + +uti_plot_show() #show all graphs (blocks script execution; close all graph windows to proceed) +print('done') diff --git a/env/python/srwpy/uti_mtrl.py b/env/python/srwpy/uti_mtrl.py new file mode 100644 index 00000000..706f40fb --- /dev/null +++ b/env/python/srwpy/uti_mtrl.py @@ -0,0 +1,413 @@ +"""Material Utilities Module + + +Modules: + + calc_refl_arr + calc_coated_refl_arr + calc_multilayer_refl_arr + add_mat + + +.. moduleauthor:: Nathan Whittington +.. moduleauthor:: Roman Chernikov +""" +# updated 22-06-2026 + +from __future__ import print_function # Python 2.7 compatibility + +import inspect +from array import array + +import numpy as np + + +xraydb = None +xraydb_found = False +xraydb_version = '' + +try: + import xraydb + xraydb_found = True + xraydb_version = getattr(xraydb, '__version__', '') +except Exception as exc: + print("Warning: 'xraydb' can not be loaded: {}. Reflectivity is set to 1.".format(exc)) + + +def _xraydb_function(_name): + if not xraydb_found: + return None + return getattr(xraydb, _name, None) + + +mirror_reflectivity = _xraydb_function('mirror_reflectivity') +coated_reflectivity = _xraydb_function('coated_reflectivity') +multilayer_reflectivity = _xraydb_function('multilayer_reflectivity') +get_material = _xraydb_function('get_material') +add_material = _xraydb_function('add_material') + + +def _warn_refl_fallback(_message): + print("Warning: {} Reflectivity is set to 1.".format(_message)) + return 1 + + +def _check_xraydb_function(_func, _func_name): + if not xraydb_found: + return _warn_refl_fallback("'xraydb' is not available.") + if _func is None: + version_text = " version {}".format(xraydb_version) if xraydb_version else '' + return _warn_refl_fallback( + "'{}' is not available in xraydb{}.".format(_func_name, version_text) + ) + try: + if 'output' not in inspect.signature(_func).parameters: + version_text = " version {}".format(xraydb_version) if xraydb_version else '' + return _warn_refl_fallback( + "'{}' in xraydb{} does not support complex amplitude output.".format( + _func_name, version_text + ) + ) + except (TypeError, ValueError): + pass + return None + + +def _check_refl_parameters(_n_ph_en, _n_ang, _n_comp, _ph_en_start, _ph_en_fin, + _ph_en_scale_type, _ang_start, _ang_fin, _ang_scale_type): + try: + n_ph_en = int(_n_ph_en) + n_ang = int(_n_ang) + n_comp = int(_n_comp) + except (TypeError, ValueError): + return _warn_refl_fallback( + "Numbers of photon energy, grazing angle, and component points must be integers." + ) + if n_ph_en < 2: + return _warn_refl_fallback( + "At least two photon energy points are required by SRW reflectivity interpolation." + ) + if n_ang < 2: + return _warn_refl_fallback( + "At least two grazing angle points are required by SRW reflectivity interpolation." + ) + if n_comp not in (1, 2): + return _warn_refl_fallback("Number of reflectivity components must be 1 or 2.") + if _ph_en_scale_type not in ('lin', 'log'): + return _warn_refl_fallback( + "Invalid photon energy scale type '{}'; use 'lin' or 'log'.".format(_ph_en_scale_type) + ) + if _ang_scale_type not in ('lin', 'log'): + return _warn_refl_fallback( + "Invalid grazing angle scale type '{}'; use 'lin' or 'log'.".format(_ang_scale_type) + ) + if (_ph_en_scale_type == 'log') or (_ang_scale_type == 'log'): + return _warn_refl_fallback( + "Logarithmic sampling is not supported by SRW reflectivity interpolation." + ) + try: + if (_ph_en_start <= 0) or (_ph_en_fin <= 0): + return _warn_refl_fallback("Photon energy limits must be positive.") + if (_ang_start < 0) or (_ang_fin < 0): + return _warn_refl_fallback("Grazing angle limits can not be negative.") + if _ph_en_start == _ph_en_fin: + return _warn_refl_fallback("Photon energy limits must define a non-zero range.") + if _ang_start == _ang_fin: + return _warn_refl_fallback("Grazing angle limits must define a non-zero range.") + except TypeError: + return _warn_refl_fallback("Photon energy and grazing angle limits must be numbers.") + return None + + +def _sampling(_start, _fin, _n, _scale_type): + if _scale_type == 'lin': + return np.linspace(_start, _fin, _n) + print("Warning: logarithmic scale is not supported by SRW reflectivity interpolation yet.") + return np.logspace(np.log10(_start), np.log10(_fin), _n) + + +def _material_is_defined(_material, _density): + if _density is not None: + return True + if (get_material is None) or (_material is None): + return False + try: + return get_material(_material) is not None + except Exception: + return False + + +def _check_material(_material, _density, _description='material'): + if _material_is_defined(_material, _density): + return None + return _warn_refl_fallback( + "Density of {} '{}' was not found; specify the density or register the " + "material with add_mat().".format(_description, _material) + ) + + +def _refl_to_srw_array(_refl_s, _refl_p, _n_tot): + try: + refl = np.ravel(_refl_s).view(np.float64) + if _refl_p is not None: + refl = np.concatenate((refl, np.ravel(_refl_p).view(np.float64))) + except (TypeError, ValueError) as exc: + return _warn_refl_fallback( + "Calculated reflectivity can not be converted to an SRW array: {}.".format(exc) + ) + if _n_tot != len(refl): + return _warn_refl_fallback( + "Calculated reflectivity array length {} does not match expected length {}.".format( + len(refl), _n_tot + ) + ) + if not np.all(np.isfinite(refl)): + return _warn_refl_fallback("Calculated reflectivity contains invalid values.") + return array('d', refl) + + +def calc_refl_arr( + _mat, + _n_ph_en=1, + _n_ang=1, + _n_comp=1, + _ph_en_start=0, + _ph_en_fin=0, + _ph_en_scale_type='lin', + _ang_start=0, + _ang_fin=0, + _ang_scale_type='lin', + _dens=None, + _roughness=0.0 +): + """Calculate complex mirror reflectivity vs photon energy and grazing angle. + + The result is a C-aligned flat array ordered by polarization, grazing angle, + photon energy, and real / imaginary component. If xraydb or required material + data are unavailable, a warning is printed and perfect reflectivity is used. + """ + fallback = _check_xraydb_function(mirror_reflectivity, 'mirror_reflectivity') + if fallback is not None: + return fallback + fallback = _check_refl_parameters( + _n_ph_en, _n_ang, _n_comp, _ph_en_start, _ph_en_fin, + _ph_en_scale_type, _ang_start, _ang_fin, _ang_scale_type + ) + if fallback is not None: + return fallback + fallback = _check_material(_mat, _dens) + if fallback is not None: + return fallback + + n_ph_en = int(_n_ph_en) + n_ang = int(_n_ang) + n_comp = int(_n_comp) + ph_en = _sampling(_ph_en_start, _ph_en_fin, n_ph_en, _ph_en_scale_type) + ang = _sampling(_ang_start, _ang_fin, n_ang, _ang_scale_type) + + try: + refl_s = mirror_reflectivity( + _mat, ang, ph_en, density=_dens, roughness=_roughness, + polarization='s', output='amplitude' + ) + refl_p = None + if n_comp == 2: + refl_p = mirror_reflectivity( + _mat, ang, ph_en, density=_dens, roughness=_roughness, + polarization='p', output='amplitude' + ) + except Exception as exc: + return _warn_refl_fallback( + "xraydb mirror reflectivity calculation for '{}' failed: {}.".format(_mat, exc) + ) + + return _refl_to_srw_array(refl_s, refl_p, n_ph_en*n_ang*n_comp*2) + + +def calc_coated_refl_arr( + _coating, + _coating_thick, + _substr, + _n_ph_en=1, + _n_ang=1, + _n_comp=1, + _ph_en_start=0, + _ph_en_fin=0, + _ph_en_scale_type='lin', + _ang_start=0, + _ang_fin=0, + _ang_scale_type='lin', + _coating_dens=None, + _substr_dens=None, + _binder=None, + _binder_dens=None, + _binder_thick=0.0, + _surf_roughness=0.0, + _substr_roughness=0.0, +): + """Calculate complex reflectivity of a coated mirror.""" + fallback = _check_xraydb_function(coated_reflectivity, 'coated_reflectivity') + if fallback is not None: + return fallback + fallback = _check_refl_parameters( + _n_ph_en, _n_ang, _n_comp, _ph_en_start, _ph_en_fin, + _ph_en_scale_type, _ang_start, _ang_fin, _ang_scale_type + ) + if fallback is not None: + return fallback + if _coating_thick is None: + return _warn_refl_fallback("Coating thickness is not specified.") + fallback = _check_material(_coating, _coating_dens, 'coating material') + if fallback is not None: + return fallback + fallback = _check_material(_substr, _substr_dens, 'substrate material') + if fallback is not None: + return fallback + if _binder is not None: + fallback = _check_material(_binder, _binder_dens, 'binder material') + if fallback is not None: + return fallback + + n_ph_en = int(_n_ph_en) + n_ang = int(_n_ang) + n_comp = int(_n_comp) + ph_en = _sampling(_ph_en_start, _ph_en_fin, n_ph_en, _ph_en_scale_type) + ang = _sampling(_ang_start, _ang_fin, n_ang, _ang_scale_type) + kwargs = dict( + coating_dens=_coating_dens, + surface_roughness=_surf_roughness, + substrate_dens=_substr_dens, + substrate_roughness=_substr_roughness, + binder=_binder, + binder_thick=_binder_thick, + binder_dens=_binder_dens, + output='amplitude', + ) + + try: + refl_s = coated_reflectivity( + _coating, _coating_thick, _substr, ang, ph_en, + polarization='s', **kwargs + ) + refl_p = None + if n_comp == 2: + refl_p = coated_reflectivity( + _coating, _coating_thick, _substr, ang, ph_en, + polarization='p', **kwargs + ) + except Exception as exc: + return _warn_refl_fallback( + "xraydb coated reflectivity calculation failed: {}.".format(exc) + ) + + return _refl_to_srw_array(refl_s, refl_p, n_ph_en*n_ang*n_comp*2) + + +def calc_multilayer_refl_arr( + _stackup, + _thickness, + _substr, + _n_periods, + _n_ph_en=1, + _n_ang=1, + _n_comp=1, + _ph_en_start=0, + _ph_en_fin=0, + _ph_en_scale_type='lin', + _ang_start=0, + _ang_fin=0, + _ang_scale_type='lin', + _dens=None, + _substr_dens=None, + _substr_rough=0, + _surf_rough=0, +): + """Calculate complex reflectivity of a periodic multilayer mirror.""" + fallback = _check_xraydb_function(multilayer_reflectivity, 'multilayer_reflectivity') + if fallback is not None: + return fallback + fallback = _check_refl_parameters( + _n_ph_en, _n_ang, _n_comp, _ph_en_start, _ph_en_fin, + _ph_en_scale_type, _ang_start, _ang_fin, _ang_scale_type + ) + if fallback is not None: + return fallback + if (_stackup is None) or (_thickness is None): + return _warn_refl_fallback("Multilayer materials and thicknesses must be specified.") + try: + n_layers = len(_stackup) + if n_layers != len(_thickness): + return _warn_refl_fallback( + "Number of multilayer materials does not match number of thicknesses." + ) + except TypeError: + return _warn_refl_fallback("Multilayer materials and thicknesses must be sequences.") + try: + n_periods = int(_n_periods) + except (TypeError, ValueError): + return _warn_refl_fallback("Number of multilayer periods must be an integer.") + if n_periods < 1: + return _warn_refl_fallback("Number of multilayer periods must be positive.") + if _dens is not None: + try: + if len(_dens) != n_layers: + return _warn_refl_fallback( + "Number of multilayer densities does not match number of materials." + ) + except TypeError: + return _warn_refl_fallback("Multilayer densities must be a sequence.") + + densities = [None]*len(_stackup) if _dens is None else _dens + for material, density in zip(_stackup, densities): + fallback = _check_material(material, density, 'multilayer material') + if fallback is not None: + return fallback + fallback = _check_material(_substr, _substr_dens, 'substrate material') + if fallback is not None: + return fallback + + n_ph_en = int(_n_ph_en) + n_ang = int(_n_ang) + n_comp = int(_n_comp) + ph_en = _sampling(_ph_en_start, _ph_en_fin, n_ph_en, _ph_en_scale_type) + ang = _sampling(_ang_start, _ang_fin, n_ang, _ang_scale_type) + kwargs = dict( + n_periods=n_periods, + density=_dens, + substrate_density=_substr_dens, + substrate_rough=_substr_rough, + surface_rough=_surf_rough, + output='amplitude', + ) + + try: + refl_s = multilayer_reflectivity( + _stackup, _thickness, _substr, ang, ph_en, + polarization='s', **kwargs + ) + refl_p = None + if n_comp == 2: + refl_p = multilayer_reflectivity( + _stackup, _thickness, _substr, ang, ph_en, + polarization='p', **kwargs + ) + except Exception as exc: + return _warn_refl_fallback( + "xraydb multilayer reflectivity calculation failed: {}.".format(exc) + ) + + return _refl_to_srw_array(refl_s, refl_p, n_ph_en*n_ang*n_comp*2) + + +def add_mat(name, formula, density, categories=None): + """Add a material to the user-local xraydb material database.""" + if add_material is None: + print("Warning: 'xraydb.add_material' is not available. Material was not added.") + return + if density is None: + print("Warning: density of '{}' is not specified. Material was not added.".format(formula)) + return + try: + add_material(name, formula, density, categories) + except Exception as exc: + print("Warning: material '{}' was not added: {}".format(name, exc)) From 60f4c48e258e082cc4a818bc9591d0144dd0de6e Mon Sep 17 00:00:00 2001 From: Roman Chernikov Date: Mon, 22 Jun 2026 22:23:47 -0400 Subject: [PATCH 2/6] Refractive properties calculation with xraydb --- env/python/srwpy/uti_mtrl.py | 81 +++++++++++++++++++++++++++++++++++- 1 file changed, 80 insertions(+), 1 deletion(-) diff --git a/env/python/srwpy/uti_mtrl.py b/env/python/srwpy/uti_mtrl.py index 706f40fb..ab1a4d72 100644 --- a/env/python/srwpy/uti_mtrl.py +++ b/env/python/srwpy/uti_mtrl.py @@ -6,6 +6,7 @@ calc_refl_arr calc_coated_refl_arr calc_multilayer_refl_arr + calc_delta_atten_len add_mat @@ -31,7 +32,7 @@ xraydb_found = True xraydb_version = getattr(xraydb, '__version__', '') except Exception as exc: - print("Warning: 'xraydb' can not be loaded: {}. Reflectivity is set to 1.".format(exc)) + print("Warning: 'xraydb' can not be loaded: {}. Reflectivity is set to 1 and refractive optical elements use vacuum properties.".format(exc)) def _xraydb_function(_name): @@ -43,6 +44,7 @@ def _xraydb_function(_name): mirror_reflectivity = _xraydb_function('mirror_reflectivity') coated_reflectivity = _xraydb_function('coated_reflectivity') multilayer_reflectivity = _xraydb_function('multilayer_reflectivity') +xray_delta_beta = _xraydb_function('xray_delta_beta') get_material = _xraydb_function('get_material') add_material = _xraydb_function('add_material') @@ -146,6 +148,34 @@ def _check_material(_material, _density, _description='material'): ) +def _refractive_fallback(_ph_en, _message): + print("Warning: {} Refractive index decrement is set to 0 and attenuation is disabled.".format(_message)) + try: + energy = np.asarray(_ph_en, dtype=float) + except (TypeError, ValueError): + return 0.0, 1.e+23 + if energy.ndim == 0: + return 0.0, 1.e+23 + return array('d', [0.0]*energy.size), array('d', [1.e+23]*energy.size) + + +def _resolve_material(_material, _density): + if _density is not None: + try: + density = float(_density) + except (TypeError, ValueError): + return None + if (not np.isfinite(density)) or (density <= 0): + return None + return _material, density + if get_material is None: + return None + try: + return get_material(_material) + except Exception: + return None + + def _refl_to_srw_array(_refl_s, _refl_p, _n_tot): try: refl = np.ravel(_refl_s).view(np.float64) @@ -399,6 +429,55 @@ def calc_multilayer_refl_arr( return _refl_to_srw_array(refl_s, refl_p, n_ph_en*n_ang*n_comp*2) +def calc_delta_atten_len(_mat, _ph_en, _dens=None): + """Calculate material properties used by SRW refractive optical elements. + + :param _mat: material name or chemical formula + :param _ph_en: photon energy [eV], or a sequence of photon energies + :param _dens: material density [g/cm^3]; if omitted, use the xraydb material database + :return: refractive index decrement and attenuation length [m] + + Scalar photon energy returns two floats. A sequence returns two ``array('d')`` + objects suitable for spectrally-dependent SRW transmission elements and CRLs. + """ + if (not xraydb_found) or (xray_delta_beta is None): + return _refractive_fallback(_ph_en, "'xraydb.xray_delta_beta' is unavailable.") + + try: + energy = np.asarray(_ph_en, dtype=float) + except (TypeError, ValueError): + return _refractive_fallback(_ph_en, "Photon energy must be numeric.") + if (energy.size < 1) or (not np.all(np.isfinite(energy))) or np.any(energy <= 0): + return _refractive_fallback(_ph_en, "Photon energy values must be finite and positive.") + + material = _resolve_material(_mat, _dens) + if material is None: + return _refractive_fallback( + _ph_en, + "Density of material '{}' was not found; specify _dens or register the material with add_mat().".format(_mat) + ) + formula, density = material + + try: + energy_arg = float(energy) if energy.ndim == 0 else energy + delta, _beta, atten_len_cm = xray_delta_beta(formula, density, energy_arg) + delta = np.asarray(delta, dtype=float) + atten_len_m = 0.01*np.asarray(atten_len_cm, dtype=float) + except Exception as exc: + return _refractive_fallback( + _ph_en, "xraydb refractive-property calculation for '{}' failed: {}.".format(_mat, exc) + ) + + if (not np.all(np.isfinite(delta))) or (not np.all(np.isfinite(atten_len_m))) or np.any(atten_len_m <= 0): + return _refractive_fallback( + _ph_en, "Calculated refractive properties for '{}' contain invalid values.".format(_mat) + ) + + if energy.ndim == 0: + return float(delta), float(atten_len_m) + return array('d', delta.ravel()), array('d', atten_len_m.ravel()) + + def add_mat(name, formula, density, categories=None): """Add a material to the user-local xraydb material database.""" if add_material is None: From 28db12ad044b0ac3e24ad39dcc332b3931837924 Mon Sep 17 00:00:00 2001 From: Roman Chernikov Date: Mon, 22 Jun 2026 22:24:25 -0400 Subject: [PATCH 3/6] Add example 22 with material properties --- .../examples/SRWLIB_Example09-materials.py | 346 ------------------ env/python/srwpy/examples/SRWLIB_Example22.py | 242 ++++++++++++ .../srwpy/examples/SRWLIB_ExamplesRunAll.py | 1 + 3 files changed, 243 insertions(+), 346 deletions(-) delete mode 100644 env/python/srwpy/examples/SRWLIB_Example09-materials.py create mode 100644 env/python/srwpy/examples/SRWLIB_Example22.py diff --git a/env/python/srwpy/examples/SRWLIB_Example09-materials.py b/env/python/srwpy/examples/SRWLIB_Example09-materials.py deleted file mode 100644 index 894afccb..00000000 --- a/env/python/srwpy/examples/SRWLIB_Example09-materials.py +++ /dev/null @@ -1,346 +0,0 @@ -############################################################################# -# SRWLIB Example#9_Materials: Simulating propagation of a Gaussian X-ray beam -# through a Beamline Containing Mirrors made of different materials -# v 0.04 -############################################################################# - -from __future__ import print_function #Python 2.7 compatibility - -try: #OC15112022 - import sys - sys.path.append('../') - from srwlib import * - from uti_plot import * - from uti_mtrl import * -except: - from srwpy.srwlib import * - from srwpy.uti_plot import * - from srwpy.uti_mtrl import * -import os -import copy - -print('SRWLIB Python Example # 9 with material properties:') -print('Simulating propagation of a coherent Gaussian X-ray beam through Pt and B4C mirrors') - - - -#******************* Reflectivity Parameters - -_n_ph_en=100 -_n_ang=1000 -_n_comp=2 -_ph_en_start=12000 -_ph_en_fin=13000 -_ph_en_scale_type="lin" -_ang_start=0.001 -_ang_fin=0.01 -_ang_scale_type="lin" -avgEn = 12400 - - -Pt_refl=calc_refl_arr('Pt', _n_ang=_n_ang, _n_ph_en=_n_ph_en, _n_comp=_n_comp, _ph_en_start=_ph_en_start, _ph_en_fin=_ph_en_fin, _ph_en_scale_type=_ph_en_scale_type, - _ang_start=_ang_start, _ang_fin=_ang_fin, _ang_scale_type=_ang_scale_type) -B4C_refl=calc_refl_arr('B4C', _dens=2.5, _n_ang=_n_ang, _n_ph_en=_n_ph_en, _n_comp=_n_comp, _ph_en_start=_ph_en_start, _ph_en_fin=_ph_en_fin, _ph_en_scale_type=_ph_en_scale_type, - _ang_start=_ang_start, _ang_fin=_ang_fin, _ang_scale_type=_ang_scale_type) - - - -#**********************Input Parameters and Structures -#***********Folder and Data File Names -strDataFolderName = 'data_example_09-materials' #data sub-folder name -strIntOutFileName = 'ex09_res_int_in.dat' #initial wavefront intensity distribution output file name -strPhOutFileName = 'ex09_res_phase_in.dat' #initial wavefront phase output file name -strIntPropOutFileName = 'ex09_res_int_prop.dat' #propagated wavefront intensity distribution output file name -strPhasePropOutFileName = 'ex09_res_phase_prop.dat' #propagated wavefront phase output file name -strIntPropOutFileName2 = 'ex09_res_int_prop_Pt.dat' #propagated wavefront intensity distribution output file name -strPhPropOutFileName2 = 'ex09_res_phase_prop_Pt.dat' #propagated wavefront phase output file name -strIntPropOutFileName3 = 'ex09_res_int_prop_B4C.dat' #propagated wavefront intensity distribution output file name -strPhPropOutFileName3 = 'ex09_res_phase_prop_B4C.dat' #propagated wavefront phase output file name - - -#***********Gaussian Beam Source -GsnBm = SRWLGsnBm() #Gaussian Beam structure (just parameters) -GsnBm.x = 0 #Transverse Positions of Gaussian Beam Center at Waist [m] -GsnBm.y = 0 -GsnBm.z = 0 #Longitudinal Position of Waist [m] -GsnBm.xp = 0 #Average Angles of Gaussian Beam at Waist [rad] -GsnBm.yp = 0 -GsnBm.avgPhotEn = avgEn #Photon Energy [eV] -GsnBm.pulseEn = 0.001 #Energy per Pulse [J] - to be corrected -GsnBm.repRate = 1 #Rep. Rate [Hz] - to be corrected -GsnBm.polar = 1 #1- linear horizontal -GsnBm.sigX = 23e-06/2.35 #Horiz. RMS size at Waist [m] -GsnBm.sigY = GsnBm.sigX #Vert. RMS size at Waist [m] - -constConvRad = 1.23984186e-06/(4*3.1415926536) -rmsAngDiv = constConvRad/(GsnBm.avgPhotEn*GsnBm.sigX) #RMS angular divergence [rad] -print('RMS Source Size:', round(GsnBm.sigX*1.e+06, 3), 'micro-m; RMS Divergence:', round(rmsAngDiv*1.e+06, 3), 'micro-rad') - -GsnBm.sigT = 10e-15 #Pulse duration [fs] (not used?) -GsnBm.mx = 0 #Transverse Gauss-Hermite Mode Orders -GsnBm.my = 0 - -#***********Initial Wavefront -wfr = SRWLWfr() #Initial Electric Field Wavefront -wfr.allocate(1, 100, 100) #Numbers of points vs Photon Energy (1), Horizontal and Vertical Positions (dummy) -wfr.mesh.zStart = 25 #Longitudinal Position [m] at which initial Electric Field has to be calculated, i.e. the position of the first optical element -wfr.mesh.eStart = GsnBm.avgPhotEn #Initial Photon Energy [eV] -wfr.mesh.eFin = GsnBm.avgPhotEn #Final Photon Energy [eV] - -wfr.unitElFld = 2 #Electric field units: 0- arbitrary, 1- sqrt(Phot/s/0.1%bw/mm^2), 2- sqrt(J/eV/mm^2) or sqrt(W/mm^2), depending on representation (freq. or time) - -distSrc_VFM = wfr.mesh.zStart - GsnBm.z -#Horizontal and Vertical Position Range for the Initial Wavefront calculation -#can be used to simulate the First Aperture -firstHorAp = 8.*rmsAngDiv*distSrc_VFM #[m] -firstVertAp = firstHorAp #[m] - -wfr.mesh.xStart = -0.5*firstHorAp #Initial Horizontal Position [m] -wfr.mesh.xFin = 0.5*firstHorAp #Final Horizontal Position [m] -wfr.mesh.yStart = -0.5*firstVertAp #Initial Vertical Position [m] -wfr.mesh.yFin = 0.5*firstVertAp #Final Vertical Position [m]s - -sampFactNxNyForProp = 4 #sampling factor for adjusting nx, ny (effective if > 0) -arPrecPar = [sampFactNxNyForProp] - -wfr.partBeam.partStatMom1.x = GsnBm.x #Some information about the source in the Wavefront structure -wfr.partBeam.partStatMom1.y = GsnBm.y -wfr.partBeam.partStatMom1.z = GsnBm.z -wfr.partBeam.partStatMom1.xp = GsnBm.xp -wfr.partBeam.partStatMom1.yp = GsnBm.yp - - -#**********************Calculating Initial Wavefront and extracting Intensity: -srwl.CalcElecFieldGaussian(wfr, GsnBm, arPrecPar) -arI0 = array('f', [0]*wfr.mesh.nx*wfr.mesh.ny) #"flat" array to take 2D intensity data -srwl.CalcIntFromElecField(arI0, wfr, 6, 0, 3, wfr.mesh.eStart, 0, 0) #extracts intensity -srwl_uti_save_intens_ascii(arI0, wfr.mesh, os.path.join(os.getcwd(), strDataFolderName, strIntOutFileName), 0) - -arI0x = array('f', [0]*wfr.mesh.nx) #array to take 1D intensity data -srwl.CalcIntFromElecField(arI0x, wfr, 6, 0, 1, wfr.mesh.eStart, 0, 0) #extracts intensity - -arP0 = array('d', [0]*wfr.mesh.nx*wfr.mesh.ny) #"flat" array to take 2D phase data (note it should be 'd') -srwl.CalcIntFromElecField(arP0, wfr, 0, 4, 3, wfr.mesh.eStart, 0, 0) #extracts radiation phase -srwl_uti_save_intens_ascii(arP0, wfr.mesh, os.path.join(os.getcwd(), strDataFolderName, strPhOutFileName), 0, ['', 'Horizontal Position', 'Vertical Position', 'Phase'], _arUnits=['', 'm', 'm', 'rad']) - -mesh0 = deepcopy(wfr.mesh) - - - -#***********Optical Elements and Propagation Parameters -#Sequence of Optical Elements: -# -# -# -# -# -# - - - -distVFM_HFM = 1 #Distance from VFM to HFM [m] -distHFM_Samp = 1 #Distance from HFM to Sample [m] - - -lenKB = 0.1 #Length of VFM and HFM (same for each) [m] -angKB = 0.004 #b4c_ca #2.521e-3 #Incident Angle of VFM and HFM [rad] - - - -#Aperture of KB system -opApKB = SRWLOptA('r', 'a', lenKB*angKB, lenKB*angKB) - - -#VFM simulated by Ideal Lens: -#opVFM = SRWLOptL(_Fy=(distSrc_M1 + distM1_VFM)*(distVFM_HFM + distHFM_Samp)/(distSrc_M1 + distM1_VFM + distVFM_HFM + distHFM_Samp)) -#VFM simulated by Extended Elliptical Mirror: -opVFM = SRWLOptMirEl(_p=(distSrc_VFM), _q=(distVFM_HFM + distHFM_Samp), _ang_graz=angKB, _size_tang=lenKB, _size_sag=10.e-03, - _nvx=0, _nvy=cos(angKB), _nvz=-sin(angKB), _tvx=0, _tvy=-sin(angKB)) - -#Drift from VFM to HFM -opDrVFM_HFM = SRWLOptD(distVFM_HFM) - -#VFM simulated by Extended Elliptical Mirror: -opHFM = SRWLOptMirEl(_p=(distSrc_VFM + distVFM_HFM), _q=distHFM_Samp, _ang_graz=angKB, _size_tang=lenKB, _size_sag=10.e-03, - _nvx=cos(angKB), _nvy=0, _nvz=-sin(angKB), _tvx=-sin(angKB), _tvy=0) - -#Drift from HFM to Sample -opDrHFM_Samp = SRWLOptD(distHFM_Samp) - - -#Wavefront Propagation Parameters: -#[0]: Auto-Resize (1) or not (0) Before propagation -#[1]: Auto-Resize (1) or not (0) After propagation -#[2]: Relative Precision for propagation with Auto-Resizing (1. is nominal) -#[3]: Allow (1) or not (0) for semi-analytical treatment of the quadratic (leading) phase terms at the propagation -#[4]: Do any Resizing on Fourier side, using FFT, (1) or not (0) -#[5]: Horizontal Range modification factor at Resizing (1. means no modification) -#[6]: Horizontal Resolution modification factor at Resizing -#[7]: Vertical Range modification factor at Resizing -#[8]: Vertical Resolution modification factor at Resizing -#[9]: Type of wavefront Shift before Resizing (not yet implemented) -#[10]: New Horizontal wavefront Center position after Shift (not yet implemented) -#[11]: New Vertical wavefront Center position after Shift (not yet implemented) -# [0][1][2] [3][4] [5] [6] [7] [8] [9][10][11] -#prParInit = [0, 0, 1., 1, 0, 2., 5., 2., 3., 0, 0, 0] -prParInit = [0, 0, 1., 1, 0, 2., 2., 2., 2., 0, 0, 0] -prPar0 = [0, 0, 1., 1, 0, 1., 2, 1., 2, 0, 0, 0] -prPar1 = [0, 0, 1., 1, 0, 1, 2, 0.1, 2, 0, 0, 0] -#prParPost = [0, 0, 1., 1, 0, 0.06, 3., 0.1, 2., 0, 0, 0] - -#NOTE: in this case of simulation, it can be enough to define the precision parameters only Before and After -#the propagation through the entire Beamline. However, if necessary, different propagation parameters can be specified -#for each optical element. -#The optimal values of propagation parameters may depend on photon energy and optical layout. - -#"Beamline" - a sequenced Container of Optical Elements (together with the corresponding wavefront propagation parameters, -#and the "post-propagation" wavefront resizing parameters for better viewing). -# optBL = SRWLOptC([opVFM, opTrErVFM, opDrVFM_HFM, opHFM, opTrErHFM, opDrHFM_Samp], -# [prPar0, prPar0, prPar0, prPar0, prPar0, prParPost]) -optBL = SRWLOptC([opApKB, opVFM, opDrVFM_HFM, opHFM, opDrHFM_Samp], - [prParInit, prPar0, prPar0, prPar0, prPar1]) - - -# # Beamline with material properties -# #VFM simulated by Ideal Lens: -# #opVFM = SRWLOptL(_Fy=(distSrc_M1 + distM1_VFM)*(distVFM_HFM + distHFM_Samp)/(distSrc_M1 + distM1_VFM + distVFM_HFM + distHFM_Samp)) -# #VFM simulated by Extended Elliptical Mirror: -opVFM_Pt = SRWLOptMirEl(_p=(distSrc_VFM), _q=(distVFM_HFM + distHFM_Samp), _ang_graz=angKB, _size_tang=lenKB, _size_sag=10.e-03, - _nvx=0, _nvy=cos(angKB), _nvz=-sin(angKB), _tvx=0, _tvy=-sin(angKB), - _refl=Pt_refl, _n_ang=_n_ang, _n_ph_en=_n_ph_en, _n_comp=_n_comp, _ph_en_start=_ph_en_start, _ph_en_fin=_ph_en_fin, _ph_en_scale_type=_ph_en_scale_type, - _ang_start=_ang_start, _ang_fin=_ang_fin, _ang_scale_type=_ang_scale_type) -# # opTrErVFM_Pt = srwl_opt_setup_surf_height_1d(heightProfData, _dim='y', _ang=angKB, _amp_coef=1) - -# #HFM simulated by Ideal Lens: -# #opHFM = SRWLOptL(_Fx=(distSrc_M1 + distM1_VFM + distVFM_HFM)*distHFM_Samp/(distSrc_M1 + distM1_VFM + distVFM_HFM + distHFM_Samp)) -# #VFM simulated by Extended Elliptical Mirror: -opHFM_Pt = SRWLOptMirEl(_p=(distSrc_VFM + distVFM_HFM), _q=distHFM_Samp, _ang_graz=angKB, _size_tang=lenKB, _size_sag=10.e-03, - _nvx=cos(angKB), _nvy=0, _nvz=-sin(angKB), _tvx=-sin(angKB), _tvy=0, - _refl=Pt_refl, _n_ang=_n_ang, _n_ph_en=_n_ph_en, _n_comp=_n_comp, _ph_en_start=_ph_en_start, _ph_en_fin=_ph_en_fin, _ph_en_scale_type=_ph_en_scale_type, - _ang_start=_ang_start, _ang_fin=_ang_fin, _ang_scale_type=_ang_scale_type) -# opTrErHFM_Pt = srwl_opt_setup_surf_height_1d(heightProfData, _dim='x', _ang=angKB, _amp_coef=1) - -# #"Beamline" - a sequenced Container of Optical Elements (together with the corresponding wavefront propagation parameters, -# #and the "post-propagation" wavefront resizing parameters for better viewing). -# # optBL_Pt = SRWLOptC([opVFM_Pt, opTrErVFM_Pt, opDrVFM_HFM, opHFM_Pt, opTrErHFM_Pt, opDrHFM_Samp], -# # [prPar0, prPar0, prPar0, prPar0, prPar0, prParPost]) -optBL_Pt = SRWLOptC([opApKB, opVFM_Pt, opDrVFM_HFM, opHFM_Pt, opDrHFM_Samp], - [prParInit, prPar0, prPar0, prPar0, prPar1]) - - - -# Boron Carbide Mirrors -opVFM_B4C = SRWLOptMirEl(_p=(distSrc_VFM), _q=(distVFM_HFM + distHFM_Samp), _ang_graz=angKB, _size_tang=lenKB, _size_sag=10.e-03, - _nvx=0, _nvy=cos(angKB), _nvz=-sin(angKB), _tvx=0, _tvy=-sin(angKB), - _refl=B4C_refl, _n_ang=_n_ang, _n_ph_en=_n_ph_en, _n_comp=_n_comp, _ph_en_start=_ph_en_start, _ph_en_fin=_ph_en_fin, _ph_en_scale_type=_ph_en_scale_type, - _ang_start=_ang_start, _ang_fin=_ang_fin, _ang_scale_type=_ang_scale_type) - -opHFM_B4C = SRWLOptMirEl(_p=(distSrc_VFM + distVFM_HFM), _q=distHFM_Samp, _ang_graz=angKB, _size_tang=lenKB, _size_sag=10.e-03, - _nvx=cos(angKB), _nvy=0, _nvz=-sin(angKB), _tvx=-sin(angKB), _tvy=0, - _refl=B4C_refl, _n_ang=_n_ang, _n_ph_en=_n_ph_en, _n_comp=_n_comp, _ph_en_start=_ph_en_start, _ph_en_fin=_ph_en_fin, _ph_en_scale_type=_ph_en_scale_type, - _ang_start=_ang_start, _ang_fin=_ang_fin, _ang_scale_type=_ang_scale_type) - -optBL_B4C = SRWLOptC([opApKB, opVFM_B4C, opDrVFM_HFM, opHFM_B4C, opDrHFM_Samp], - [prParInit, prPar0, prPar0, prPar0, prPar1]) - - - - -#***********Wavefront Propagation -arI1 = None; arI1x = None; arI1y = None -arP1 = None -mesh1 = None - - -#Duplicating wavefront (by re-creating all objects/arrays): -wfr1 = copy.deepcopy(wfr) -wfr2 = copy.deepcopy(wfr) -wfr3 = copy.deepcopy(wfr) - - - -print(' Propagating Wavefront (through two elliptical mirrors)... ', end='') -srwl.PropagElecField(wfr1, optBL) -print('done') -print(' Saving resulting Intensity and Phase data to files ... ', end='') -mesh1 = deepcopy(wfr1.mesh) -arI1 = array('f', [0]*mesh1.nx*mesh1.ny) #"flat" array to take 2D intensity data -srwl.CalcIntFromElecField(arI1, wfr1, 6, 0, 3, mesh1.eStart, 0, 0) #extracts intensity -srwl_uti_save_intens_ascii(arI1, mesh1, os.path.join(os.getcwd(), strDataFolderName, strIntPropOutFileName), 0) - -arI1x = array('f', [0]*mesh1.nx) #array to take 1D intensity data -srwl.CalcIntFromElecField(arI1x, wfr1, 6, 0, 1, mesh1.eStart, 0, 0) #extracts intensity -arI1y = array('f', [0]*mesh1.ny) #array to take 1D intensity data -srwl.CalcIntFromElecField(arI1y, wfr1, 6, 0, 2, mesh1.eStart, 0, 0) #extracts intensity -arP1 = array('d', [0]*mesh1.nx*mesh1.ny) #"flat" array to take 2D phase data (note it should be 'd') -srwl.CalcIntFromElecField(arP1, wfr1, 0, 4, 3, mesh1.eStart, 0, 0) #extracts radiation phase -srwl_uti_save_intens_ascii(arP1, mesh1, os.path.join(os.getcwd(), strDataFolderName, strPhasePropOutFileName), 0, ['', 'Horizontal Position', 'Vertical Position', 'Phase'], _arUnits=['', 'm', 'm', 'rad']) -del wfr1 -print('done') - - -print(' Propagating Wavefront (through 2 Platinum Mirrors) ... ', end='') -srwl.PropagElecField(wfr2, optBL_Pt) -print('done') -print(' Saving resulting Intensity and Phase data to files ... ', end='') -mesh2 = deepcopy(wfr2.mesh) -arI2 = array('f', [0]*mesh2.nx*mesh2.ny) #"flat" array to take 2D intensity data -srwl.CalcIntFromElecField(arI2, wfr2, 6, 0, 3, mesh2.eStart, 0, 0) #extracts intensity -srwl_uti_save_intens_ascii(arI2, mesh2, os.path.join(os.getcwd(), strDataFolderName, strIntPropOutFileName2), 0) - -arI2x = array('f', [0]*mesh2.nx) #array to take 1D intensity data -srwl.CalcIntFromElecField(arI2x, wfr2, 6, 0, 1, mesh2.eStart, 0, 0) #extracts intensity -arI2y = array('f', [0]*mesh2.ny) #array to take 1D intensity data -srwl.CalcIntFromElecField(arI2y, wfr2, 6, 0, 2, mesh2.eStart, 0, 0) #extracts intensity -arP2 = array('d', [0]*mesh2.nx*mesh2.ny) #"flat" array to take 2D phase data (note it should be 'd') -srwl.CalcIntFromElecField(arP2, wfr2, 0, 4, 3, mesh2.eStart, 0, 0) #extracts radiation phase -srwl_uti_save_intens_ascii(arP2, mesh2, os.path.join(os.getcwd(), strDataFolderName, strPhPropOutFileName2), 0, ['', 'Horizontal Position', 'Vertical Position', 'Phase'], _arUnits=['', 'm', 'm', 'rad']) -del wfr2 -print('done') - -print(' Propagating Wavefront (through 2 B4C Mirrors) ... ', end='') -srwl.PropagElecField(wfr3, optBL_B4C) -print('done') -print(' Saving resulting Intensity and Phase data to files ... ', end='') -mesh3 = deepcopy(wfr3.mesh) -arI3 = array('f', [0]*mesh3.nx*mesh3.ny) -srwl.CalcIntFromElecField(arI3, wfr3, 6, 0, 3, mesh3.eStart, 0, 0) -srwl_uti_save_intens_ascii(arI3, mesh3, os.path.join(os.getcwd(), strDataFolderName, strIntPropOutFileName3), 0) - -arP3 = array('d', [0]*mesh3.nx*mesh3.ny) -srwl.CalcIntFromElecField(arP3, wfr3, 0, 4, 3, mesh3.eStart, 0, 0) -srwl_uti_save_intens_ascii(arP3, mesh3, os.path.join(os.getcwd(), strDataFolderName, strPhPropOutFileName3), 0, ['', 'Horizontal Position', 'Vertical Position', 'Phase'], _arUnits=['', 'm', 'm', 'rad']) -del wfr3 -print('done') - -#**********************Plotting results (requires 3rd party graphics package) -print(' Plotting the results (blocks script execution; close any graph windows to proceed) ... ', end='') -plotMesh0x = [1000*mesh0.xStart, 1000*mesh0.xFin, mesh0.nx] -plotMesh0y = [1000*mesh0.yStart, 1000*mesh0.yFin, mesh0.ny] -uti_plot2d1d(arI0, plotMesh0x, plotMesh0y, labels=['Horizontal Position [mm]', 'Vertical Position [mm]', 'Intensity Before Propagation']) - - -plotMesh1x = [1e+06*mesh1.xStart, 1e+06*mesh1.xFin, mesh1.nx] -plotMesh1y = [1e+06*mesh1.yStart, 1e+06*mesh1.yFin, mesh1.ny] -uti_plot2d1d(arI1, plotMesh1x, plotMesh1y, labels=['Horizontal Position [microns]', 'Vertical Position [microns]', 'Intensity After Propagation through 2 Elliptical Mirrors']) - -plotMesh2x = [1e+06*mesh2.xStart, 1e+06*mesh2.xFin, mesh2.nx] -plotMesh2y = [1e+06*mesh2.yStart, 1e+06*mesh2.yFin, mesh2.ny] -uti_plot2d1d(arI2, plotMesh2x, plotMesh2y, labels=['Horizontal Position [microns]', 'Vertical Position [microns]', 'Intensity After Propagation through 2 Platinum Elliptical Mirrors']) - -plotMesh3x = [1e+06*mesh3.xStart, 1e+06*mesh3.xFin, mesh3.nx] -plotMesh3y = [1e+06*mesh3.yStart, 1e+06*mesh3.yFin, mesh3.ny] -uti_plot2d1d(arI3, plotMesh3x, plotMesh3y, labels=['Horizontal Position [microns]', 'Vertical Position [microns]', 'Intensity After Propagation through 2 Boron Carbide Elliptical Mirrors']) - - - -# Calculate relative intensity after propagation through the material mirrors -print() -print("maximum intensity for perfect reflectivity: ", max(arI1)) -print("maximum intensity for Pt reflectivity: ", max(arI2)) -print("maximum intensity for B4C reflectivity: ", max(arI3)) -print("simulated relative intensity for Pt mirrors: ", max(arI2)/max(arI1)) -print("simulated relative intensity for B4C mirrors: ", max(arI3)/max(arI1)) - -uti_plot_show() #show all graphs (blocks script execution; close all graph windows to proceed) -print('done') diff --git a/env/python/srwpy/examples/SRWLIB_Example22.py b/env/python/srwpy/examples/SRWLIB_Example22.py new file mode 100644 index 00000000..16414c83 --- /dev/null +++ b/env/python/srwpy/examples/SRWLIB_Example22.py @@ -0,0 +1,242 @@ +############################################################################# +# SRWLIB Example # 22: Material-dependent reflective and refractive optics +# Simulating a Gaussian X-ray beam passing through an Aluminum filter, +# a B4C mirror, a 200 nm Pt-coated Si mirror, and an array of diamond CRLs +# v 0.01 +############################################################################# + +from __future__ import print_function + +try: #OC15112022 + import sys + sys.path.append('../') + from srwlib import * + from uti_plot import * + from uti_mtrl import * +except: + from srwpy.srwlib import * + from srwpy.uti_plot import * + from srwpy.uti_mtrl import * + +import math +import os +from copy import deepcopy + + +print('SRWLIB Python Example # 22:') +print('Simulating material-dependent reflective and refractive X-ray optics') + + +#******************* Photon Energy and Material Parameters +photon_energy = 8000. # [eV] + +n_ph_en = 101 +n_ang = 401 +n_comp = 2 +ph_en_start = 7500. +ph_en_fin = 8500. +ang_start = 0.001 +ang_fin = 0.005 +ang_graz = 0.003 + +# Aluminum filter +filter_thick = 50.e-06 # [m] +al_delta, al_atten_len = calc_delta_atten_len('Al', photon_energy) + +# B4C mirror; B4C is a formula not registered with a default xraydb density +b4c_density = 2.5 # [g/cm^3] +b4c_refl = calc_refl_arr( + 'B4C', _dens=b4c_density, + _n_ph_en=n_ph_en, _n_ang=n_ang, _n_comp=n_comp, + _ph_en_start=ph_en_start, _ph_en_fin=ph_en_fin, + _ang_start=ang_start, _ang_fin=ang_fin, +) + +# 200 nm Pt coating on a Si substrate; xraydb expects thickness in Angstroms +pt_coating_thick = 2000. # [Angstrom] +pt_si_refl = calc_coated_refl_arr( + 'Pt', pt_coating_thick, 'Si', + _n_ph_en=n_ph_en, _n_ang=n_ang, _n_comp=n_comp, + _ph_en_start=ph_en_start, _ph_en_fin=ph_en_fin, + _ang_start=ang_start, _ang_fin=ang_fin, +) + +# Diamond CRL material: elemental C at diamond density +diamond_density = 3.52 # [g/cm^3] +diamond_delta, diamond_atten_len = calc_delta_atten_len( + 'C', photon_energy, _dens=diamond_density +) + +print('Al filter: delta = {:.6g}, attenuation length = {:.6g} m'.format( + al_delta, al_atten_len +)) +print('Diamond: delta = {:.6g}, attenuation length = {:.6g} m'.format( + diamond_delta, diamond_atten_len +)) + + +#********************** Output Files +strDataFolderName = 'data_example_22' +strIntInFileName = 'ex22_res_int_in.dat' +strIntOutFileName = 'ex22_res_int_out.dat' + +if not os.path.exists(strDataFolderName): + os.makedirs(strDataFolderName) + + +#********************** Gaussian Beam Source +gsnBm = SRWLGsnBm() +gsnBm.x = 0 +gsnBm.y = 0 +gsnBm.z = 0 +gsnBm.xp = 0 +gsnBm.yp = 0 +gsnBm.avgPhotEn = photon_energy +gsnBm.pulseEn = 0.001 +gsnBm.repRate = 1 +gsnBm.polar = 1 # linear horizontal polarization +gsnBm.sigX = 20.e-06 +gsnBm.sigY = 20.e-06 +gsnBm.sigT = 10.e-15 +gsnBm.mx = 0 +gsnBm.my = 0 + +wfr = SRWLWfr() +wfr.allocate(1, 151, 151) +wfr.mesh.zStart = 10. +wfr.mesh.eStart = photon_energy +wfr.mesh.eFin = photon_energy +wfr.mesh.xStart = -0.3e-03 +wfr.mesh.xFin = 0.3e-03 +wfr.mesh.yStart = -0.3e-03 +wfr.mesh.yFin = 0.3e-03 +wfr.unitElFld = 2 +wfr.partBeam.partStatMom1.x = gsnBm.x +wfr.partBeam.partStatMom1.y = gsnBm.y +wfr.partBeam.partStatMom1.z = gsnBm.z +wfr.partBeam.partStatMom1.xp = gsnBm.xp +wfr.partBeam.partStatMom1.yp = gsnBm.yp + +srwl.CalcElecFieldGaussian(wfr, gsnBm, [2]) + +meshIn = deepcopy(wfr.mesh) +arIIn = array('f', [0]*meshIn.nx*meshIn.ny) +srwl.CalcIntFromElecField(arIIn, wfr, 6, 0, 3, meshIn.eStart, 0, 0) +srwl_uti_save_intens_ascii( + arIIn, meshIn, os.path.join(strDataFolderName, strIntInFileName), 0 +) + + +#********************** Optical Elements +# Uniform Aluminum filter. Amplitude transmission is exp(-thickness/(2*L)). +filter_amp = math.exp(-0.5*filter_thick/al_atten_len) +filter_opd = -al_delta*filter_thick +filter_ar_tr = array('d', [filter_amp, filter_opd]*4) +opFilter = SRWLOptT( + _nx=2, _ny=2, _rx=2.e-03, _ry=2.e-03, + _arTr=filter_ar_tr, _extTr=1 +) + +mirror_len = 0.2 +mirror_width = 5.e-03 + +# B4C plane mirror, deflecting vertically +opMirB4C = SRWLOptMirPl( + _size_tang=mirror_len, _size_sag=mirror_width, + _nvx=0, _nvy=cos(ang_graz), _nvz=-sin(ang_graz), + _tvx=0, _tvy=-sin(ang_graz), + _refl=b4c_refl, + _n_ph_en=n_ph_en, _n_ang=n_ang, _n_comp=n_comp, + _ph_en_start=ph_en_start, _ph_en_fin=ph_en_fin, + _ang_start=ang_start, _ang_fin=ang_fin, +) + +# Pt-coated Si plane mirror, deflecting horizontally +opMirPtSi = SRWLOptMirPl( + _size_tang=mirror_len, _size_sag=mirror_width, + _nvx=cos(ang_graz), _nvy=0, _nvz=-sin(ang_graz), + _tvx=-sin(ang_graz), _tvy=0, + _refl=pt_si_refl, + _n_ph_en=n_ph_en, _n_ang=n_ang, _n_comp=n_comp, + _ph_en_start=ph_en_start, _ph_en_fin=ph_en_fin, + _ang_start=ang_start, _ang_fin=ang_fin, +) + +# Array of ten 2D parabolic diamond CRLs +crl_apert = 1.e-03 +crl_r_min = 100.e-06 +crl_number = 10 +crl_wall_thick = 20.e-06 +if diamond_delta > 0: + opCRL = srwl_opt_setup_CRL( + 3, diamond_delta, diamond_atten_len, 1, + crl_apert, crl_apert, crl_r_min, crl_number, crl_wall_thick, + 0, 0, None, 0, 0, 201, 201 + ) + crl_focal_len = opCRL.Fx +else: + print('Warning: diamond refractive properties are unavailable; the CRL is disabled.') + opCRL = SRWLOptT( + _nx=2, _ny=2, _rx=crl_apert, _ry=crl_apert, + _extTr=1, _alloc_base=[1, 0] + ) + crl_focal_len = 1. + +opDrFilter_M1 = SRWLOptD(1.) +opDrM1_M2 = SRWLOptD(1.) +opDrM2_CRL = SRWLOptD(1.) +opDrCRL_Obs = SRWLOptD(crl_focal_len) + + +#********************** Propagation +ppOpt = [0, 0, 1., 1, 0, 1., 1., 1., 1., 0, 0, 0] +ppDrift = [0, 0, 1., 1, 0, 1., 1., 1., 1., 0, 0, 0] +ppCRL = [0, 0, 1., 1, 0, 1.1, 2., 1.1, 2., 0, 0, 0] +ppFinal = [0, 0, 1., 1, 0, 0.2, 2., 0.2, 2., 0, 0, 0] + +optBL = SRWLOptC( + [opFilter, opDrFilter_M1, opMirB4C, opDrM1_M2, + opMirPtSi, opDrM2_CRL, opCRL, opDrCRL_Obs], + [ppOpt, ppDrift, ppOpt, ppDrift, ppOpt, ppDrift, ppCRL, ppFinal] +) + +print(' Propagating wavefront through filter, mirrors, and diamond CRLs ... ', end='') +srwl.PropagElecField(wfr, optBL) +print('done') + +meshOut = deepcopy(wfr.mesh) +arIOut = array('f', [0]*meshOut.nx*meshOut.ny) +srwl.CalcIntFromElecField(arIOut, wfr, 6, 0, 3, meshOut.eStart, 0, 0) +srwl_uti_save_intens_ascii( + arIOut, meshOut, os.path.join(strDataFolderName, strIntOutFileName), 0 +) + + +def integrated_intensity(arI, mesh): + dx = (mesh.xFin - mesh.xStart)/(mesh.nx - 1) if mesh.nx > 1 else 1. + dy = (mesh.yFin - mesh.yStart)/(mesh.ny - 1) if mesh.ny > 1 else 1. + return sum(arI)*dx*dy + + +intIn = integrated_intensity(arIIn, meshIn) +intOut = integrated_intensity(arIOut, meshOut) +print('CRL focal length: {:.6g} m'.format(crl_focal_len)) +print('Al filter intensity transmission: {:.6g}'.format(filter_amp*filter_amp)) +print('Integrated beamline transmission: {:.6g}'.format(intOut/intIn)) + + +#********************** Plotting +uti_plot2d1d( + arIIn, + [1.e+03*meshIn.xStart, 1.e+03*meshIn.xFin, meshIn.nx], + [1.e+03*meshIn.yStart, 1.e+03*meshIn.yFin, meshIn.ny], + labels=['Horizontal Position [mm]', 'Vertical Position [mm]', 'Input Intensity'] +) +uti_plot2d1d( + arIOut, + [1.e+06*meshOut.xStart, 1.e+06*meshOut.xFin, meshOut.nx], + [1.e+06*meshOut.yStart, 1.e+06*meshOut.yFin, meshOut.ny], + labels=['Horizontal Position [microns]', 'Vertical Position [microns]', 'Focused Intensity'] +) +uti_plot_show() +print('done') diff --git a/env/python/srwpy/examples/SRWLIB_ExamplesRunAll.py b/env/python/srwpy/examples/SRWLIB_ExamplesRunAll.py index 9774106f..200bdc55 100644 --- a/env/python/srwpy/examples/SRWLIB_ExamplesRunAll.py +++ b/env/python/srwpy/examples/SRWLIB_ExamplesRunAll.py @@ -32,6 +32,7 @@ 'SRWLIB_Example18.py', 'SRWLIB_Example19.py', 'SRWLIB_Example21.py', + 'SRWLIB_Example22.py', ] for i in range(len(exFileNames)): From 80f77b726f00a23eaf7d833d1efb67e067cf0e44 Mon Sep 17 00:00:00 2001 From: Roman Chernikov Date: Mon, 22 Jun 2026 23:13:51 -0400 Subject: [PATCH 4/6] Fix reflectivity parsing and cleanup --- cpp/src/clients/python/srwlpy.cpp | 273 +++++++++++++++++++++++++----- 1 file changed, 230 insertions(+), 43 deletions(-) diff --git a/cpp/src/clients/python/srwlpy.cpp b/cpp/src/clients/python/srwlpy.cpp index b0c12214..3f965f35 100644 --- a/cpp/src/clients/python/srwlpy.cpp +++ b/cpp/src/clients/python/srwlpy.cpp @@ -24,6 +24,7 @@ #include "pyparse.h" //OC09032019 #include #include +#include #include //OCTEST_161214 //OC18022024 (commented-out) @@ -82,6 +83,7 @@ static const char strEr_BadOptZP[] = "Incorrect Optical Zone Plate structure"; static const char strEr_BadOptWG[] = "Incorrect Optical Waveguide structure"; static const char strEr_BadOptG[] = "Incorrect Optical Grating structure"; static const char strEr_BadOptT[] = "Incorrect Optical Generic Transmission structure"; +static const char strEr_BadOptR[] = "Incorrect Reflectivity Object structure"; static const char strEr_BadOptMir[] = "Incorrect Optical Mirror structure"; static const char strEr_BadOptCryst[] = "Incorrect Optical Crystal structure"; static const char strEr_BadListIntProp[] = "Incorrect list structure defining intensity distributions to be plotted after propagation"; @@ -2019,6 +2021,96 @@ void ParseSructSRWLOptT(SRWLOptT* pOpt, PyObject* oOpt, vector* pvBuf } } +// NW11072025 +/************************************************************************** + * Parses PyObject* to SRWLOptR* + * **************************************************************************/ +void ParseSructSRWLOptR(SRWLOptR* pOpt, PyObject* oOpt, vector* pvBuf) //throw(...) +{ + if((pOpt == 0) || (oOpt == 0)) throw strEr_NoObj; + + PyObject *o_tmp = 0; + Py_ssize_t sizeRefl = 0; + o_tmp = PyObject_GetAttrString(oOpt, "arRefl"); + if(o_tmp == 0) throw strEr_BadOptR; + pOpt->arRefl = (double*)GetPyArrayBuf(o_tmp, pvBuf, &sizeRefl); + if(pOpt->arRefl == 0) { Py_DECREF(o_tmp); throw strEr_BadOptR; } + Py_DECREF(o_tmp); + + o_tmp = PyObject_GetAttrString(oOpt, "reflNumPhEn"); + if(o_tmp == 0) throw strEr_BadOptR; + if(!PyNumber_Check(o_tmp)) { Py_DECREF(o_tmp); throw strEr_BadOptR; } + pOpt->reflNumPhEn = PyLong_AsLong(o_tmp); + if(PyErr_Occurred()) { Py_DECREF(o_tmp); throw strEr_BadOptR; } + Py_DECREF(o_tmp); + + o_tmp = PyObject_GetAttrString(oOpt, "reflNumAng"); + if(o_tmp == 0) throw strEr_BadOptR; + if(!PyNumber_Check(o_tmp)) { Py_DECREF(o_tmp); throw strEr_BadOptR; } + pOpt->reflNumAng = PyLong_AsLong(o_tmp); + if(PyErr_Occurred()) { Py_DECREF(o_tmp); throw strEr_BadOptR; } + Py_DECREF(o_tmp); + + o_tmp = PyObject_GetAttrString(oOpt, "reflNumComp"); + if(o_tmp == 0) throw strEr_BadOptR; + if(!PyNumber_Check(o_tmp)) { Py_DECREF(o_tmp); throw strEr_BadOptR; } + pOpt->reflNumComp = PyLong_AsLong(o_tmp); + if(PyErr_Occurred()) { Py_DECREF(o_tmp); throw strEr_BadOptR; } + Py_DECREF(o_tmp); + + o_tmp = PyObject_GetAttrString(oOpt, "reflPhEnScaleType"); + if(o_tmp == 0) throw strEr_BadOptR; + if(!PyUnicode_Check(o_tmp) && !PyBytes_Check(o_tmp)) { Py_DECREF(o_tmp); throw strEr_BadOptR; } + CopyPyStringToC(o_tmp, pOpt->reflPhEnScaleType, 3); + Py_DECREF(o_tmp); + + o_tmp = PyObject_GetAttrString(oOpt, "reflAngScaleType"); + if(o_tmp == 0) throw strEr_BadOptR; + if(!PyUnicode_Check(o_tmp) && !PyBytes_Check(o_tmp)) { Py_DECREF(o_tmp); throw strEr_BadOptR; } + CopyPyStringToC(o_tmp, pOpt->reflAngScaleType, 3); + Py_DECREF(o_tmp); + + o_tmp = PyObject_GetAttrString(oOpt, "reflPhEnStart"); + if(o_tmp == 0) throw strEr_BadOptR; + if(!PyNumber_Check(o_tmp)) { Py_DECREF(o_tmp); throw strEr_BadOptR; } + pOpt->reflPhEnStart = PyFloat_AsDouble(o_tmp); + if(PyErr_Occurred()) { Py_DECREF(o_tmp); throw strEr_BadOptR; } + Py_DECREF(o_tmp); + + o_tmp = PyObject_GetAttrString(oOpt, "reflPhEnFin"); + if(o_tmp == 0) throw strEr_BadOptR; + if(!PyNumber_Check(o_tmp)) { Py_DECREF(o_tmp); throw strEr_BadOptR; } + pOpt->reflPhEnFin = PyFloat_AsDouble(o_tmp); + if(PyErr_Occurred()) { Py_DECREF(o_tmp); throw strEr_BadOptR; } + Py_DECREF(o_tmp); + + o_tmp = PyObject_GetAttrString(oOpt, "reflAngStart"); + if(o_tmp == 0) throw strEr_BadOptR; + if(!PyNumber_Check(o_tmp)) { Py_DECREF(o_tmp); throw strEr_BadOptR; } + pOpt->reflAngStart = PyFloat_AsDouble(o_tmp); + if(PyErr_Occurred()) { Py_DECREF(o_tmp); throw strEr_BadOptR; } + Py_DECREF(o_tmp); + + o_tmp = PyObject_GetAttrString(oOpt, "reflAngFin"); + if(o_tmp == 0) throw strEr_BadOptR; + if(!PyNumber_Check(o_tmp)) { Py_DECREF(o_tmp); throw strEr_BadOptR; } + pOpt->reflAngFin = PyFloat_AsDouble(o_tmp); + if(PyErr_Occurred()) { Py_DECREF(o_tmp); throw strEr_BadOptR; } + Py_DECREF(o_tmp); + + if((pOpt->reflNumPhEn <= 0) || (pOpt->reflNumAng <= 0) || + (pOpt->reflNumComp < 1) || (pOpt->reflNumComp > 2)) throw strEr_BadOptR; + + size_t nVal = (size_t)pOpt->reflNumPhEn; + size_t nAng = (size_t)pOpt->reflNumAng; + size_t nComp = (size_t)pOpt->reflNumComp; + size_t maxSize = numeric_limits::max(); + if((nVal > maxSize/nAng) || ((nVal *= nAng) > maxSize/nComp) || + ((nVal *= nComp) > maxSize/2) || ((nVal *= 2) > maxSize/sizeof(double))) throw strEr_BadOptR; + size_t sizeReflExp = nVal*sizeof(double); + if((sizeRefl < 0) || ((size_t)sizeRefl != sizeReflExp)) throw strEr_BadOptR; +} + /************************************************************************//** * Parses PyObject* to SRWLOptMir* ***************************************************************************/ @@ -2026,6 +2118,10 @@ void ParseSructSRWLOptMir(SRWLOptMir* pOpt, PyObject* oOpt, vector* p { if((pOpt == 0) || (oOpt == 0)) throw strEr_NoObj; + pOpt->arReflObj = 0; + pOpt->nReflObj = 0; + pOpt->arReflDist = 0; + //PyObject *o_tmp = PyObject_GetAttrString(oOpt, "arRefl"); //pOpt->arRefl = (double*)GetPyArrayBuf(o_tmp, pvBuf, 0); //if(pOpt->arRefl == 0) throw strEr_BadOptMir; @@ -2228,6 +2324,55 @@ void ParseSructSRWLOptMir(SRWLOptMir* pOpt, PyObject* oOpt, vector* p pOpt->isConvex = (char)PyLong_AsLong(o_tmp); Py_DECREF(o_tmp); } + + //NW11072025 ------------------------------ + o_tmp = PyObject_GetAttrString(oOpt, "arReflObj"); + if(o_tmp == 0) + { + PyErr_Clear(); + return; + } + if(o_tmp == Py_None) + { + Py_DECREF(o_tmp); + return; + } + if(!PyList_Check(o_tmp)) { Py_DECREF(o_tmp); throw strEr_BadOptMir; } + + int nReflObj = (int)PyList_Size(o_tmp); + if(nReflObj <= 0) { Py_DECREF(o_tmp); throw strEr_BadOptMir; } + vector vReflObj(nReflObj); + for(int i=0; inpt <= 0) || (pOpt->nps <= 0)) { Py_DECREF(o_tmp); throw strEr_BadOptMir; } + + size_t nMap = (size_t)pOpt->npt; + size_t nps = (size_t)pOpt->nps; + size_t maxSize = numeric_limits::max(); + if((nMap > maxSize/nps) || ((nMap *= nps) > maxSize/sizeof(int))) { Py_DECREF(o_tmp); throw strEr_BadOptMir; } + size_t sizeReflDistExp = nMap*sizeof(int); + if((sizeReflDist < 0) || ((size_t)sizeReflDist != sizeReflDistExp)) { Py_DECREF(o_tmp); throw strEr_BadOptMir; } + for(size_t i=0; i= nReflObj)) { Py_DECREF(o_tmp); throw strEr_BadOptMir; } + } + Py_DECREF(o_tmp); + + pOpt->arReflObj = new SRWLOptR[nReflObj]; + for(int i=0; iarReflObj[i] = vReflObj[i]; + pOpt->nReflObj = nReflObj; + pOpt->arReflDist = pReflDist; + //NW11072025 ------------------------------ } /************************************************************************//** @@ -2381,6 +2526,32 @@ void ParseSructSRWLOptMirExtHyp(SRWLOptMirHyp* pOpt, PyObject* oOpt) //throw(... Py_DECREF(o_tmp); } +/************************************************************************//** + * Deallocates a parsed mirror of a specific type. + ***************************************************************************/ +template void DeallocSRWLOptMirType(void* pMir) +{ + T *pMirTyped = (T*)pMir; + if(pMirTyped->baseMir.arReflObj != 0) delete[] pMirTyped->baseMir.arReflObj; + delete pMirTyped; +} + +/************************************************************************//** + * Deallocates a parsed mirror and its reflectivity object array. + ***************************************************************************/ +void DeallocSRWLOptMir(void* pMir, const char* sType) +{ + if((pMir == 0) || (sType == 0)) return; + + if(strcmp(sType, "mirror: plane") == 0) DeallocSRWLOptMirType(pMir); + else if(strcmp(sType, "mirror: ellipsoid") == 0) DeallocSRWLOptMirType(pMir); + else if(strcmp(sType, "mirror: paraboloid") == 0) DeallocSRWLOptMirType(pMir); + else if(strcmp(sType, "mirror: toroid") == 0) DeallocSRWLOptMirType(pMir); + else if(strcmp(sType, "mirror: sphere") == 0) DeallocSRWLOptMirType(pMir); + else if(strcmp(sType, "mirror: hyperboloid") == 0) DeallocSRWLOptMirType(pMir); +} + + /************************************************************************//** * Parses PyObject* to a SRWLOptMir* ***************************************************************************/ @@ -2401,46 +2572,54 @@ void* ParseSructSRWLOptMirAll(PyObject* oOpt, char* sPyTypeName, vectorbaseMir), oOpt, pvBuf); - } - else if(strcmp(sPyTypeName, "SRWLOptMirEl") == 0) - { - pMir = new SRWLOptMirEl(); - strcat(srwOptTypeName, "ellipsoid\0"); - ParseSructSRWLOptMir(&(((SRWLOptMirEl*)pMir)->baseMir), oOpt, pvBuf); - ParseSructSRWLOptMirExtEl((SRWLOptMirEl*)pMir, oOpt); - } - else if(strcmp(sPyTypeName, "SRWLOptMirPar") == 0) - { - pMir = new SRWLOptMirPar(); - strcat(srwOptTypeName, "paraboloid\0"); - ParseSructSRWLOptMir(&(((SRWLOptMirPar*)pMir)->baseMir), oOpt, pvBuf); - ParseSructSRWLOptMirExtPar((SRWLOptMirPar*)pMir, oOpt); - } - else if(strcmp(sPyTypeName, "SRWLOptMirTor") == 0) - { - pMir = new SRWLOptMirTor(); - strcat(srwOptTypeName, "toroid\0"); - ParseSructSRWLOptMir(&(((SRWLOptMirTor*)pMir)->baseMir), oOpt, pvBuf); - ParseSructSRWLOptMirExtTor((SRWLOptMirTor*)pMir, oOpt); - } - else if (strcmp(sPyTypeName, "SRWLOptMirSph") == 0) + try { - pMir = new SRWLOptMirSph(); - strcat(srwOptTypeName, "sphere\0"); - ParseSructSRWLOptMir(&(((SRWLOptMirSph*)pMir)->baseMir), oOpt, pvBuf); - ParseSructSRWLOptMirExtSph((SRWLOptMirSph*)pMir, oOpt); + if(strcmp(sPyTypeName, "SRWLOptMirPl") == 0) + { + pMir = new SRWLOptMirPl(); + strcat(srwOptTypeName, "plane\0"); + ParseSructSRWLOptMir(&(((SRWLOptMirPl*)pMir)->baseMir), oOpt, pvBuf); + } + else if(strcmp(sPyTypeName, "SRWLOptMirEl") == 0) + { + pMir = new SRWLOptMirEl(); + strcat(srwOptTypeName, "ellipsoid\0"); + ParseSructSRWLOptMir(&(((SRWLOptMirEl*)pMir)->baseMir), oOpt, pvBuf); + ParseSructSRWLOptMirExtEl((SRWLOptMirEl*)pMir, oOpt); + } + else if(strcmp(sPyTypeName, "SRWLOptMirPar") == 0) + { + pMir = new SRWLOptMirPar(); + strcat(srwOptTypeName, "paraboloid\0"); + ParseSructSRWLOptMir(&(((SRWLOptMirPar*)pMir)->baseMir), oOpt, pvBuf); + ParseSructSRWLOptMirExtPar((SRWLOptMirPar*)pMir, oOpt); + } + else if(strcmp(sPyTypeName, "SRWLOptMirTor") == 0) + { + pMir = new SRWLOptMirTor(); + strcat(srwOptTypeName, "toroid\0"); + ParseSructSRWLOptMir(&(((SRWLOptMirTor*)pMir)->baseMir), oOpt, pvBuf); + ParseSructSRWLOptMirExtTor((SRWLOptMirTor*)pMir, oOpt); + } + else if (strcmp(sPyTypeName, "SRWLOptMirSph") == 0) + { + pMir = new SRWLOptMirSph(); + strcat(srwOptTypeName, "sphere\0"); + ParseSructSRWLOptMir(&(((SRWLOptMirSph*)pMir)->baseMir), oOpt, pvBuf); + ParseSructSRWLOptMirExtSph((SRWLOptMirSph*)pMir, oOpt); + } + else if(strcmp(sPyTypeName, "SRWLOptMirHyp") == 0) //TW24012024 + { + pMir = new SRWLOptMirHyp(); + strcat(srwOptTypeName, "hyperboloid\0"); + ParseSructSRWLOptMir(&(((SRWLOptMirHyp*)pMir)->baseMir), oOpt, pvBuf); + ParseSructSRWLOptMirExtHyp((SRWLOptMirHyp*)pMir, oOpt); + } } - else if(strcmp(sPyTypeName, "SRWLOptMirHyp") == 0) //TW24012024 + catch(...) { - pMir = new SRWLOptMirHyp(); - strcat(srwOptTypeName, "hyperboloid\0"); - ParseSructSRWLOptMir(&(((SRWLOptMirHyp*)pMir)->baseMir), oOpt, pvBuf); - ParseSructSRWLOptMirExtHyp((SRWLOptMirHyp*)pMir, oOpt); + DeallocSRWLOptMir(pMir, srwOptTypeName); + throw; } return pMir; @@ -2452,13 +2631,10 @@ void* ParseSructSRWLOptMirAll(PyObject* oOpt, char* sPyTypeName, vector* pvBuf) //throw(...) { if((pOpt == 0) || (oOpt == 0)) throw strEr_NoObj; + pOpt->mirSub = 0; + pOpt->mirSubType[0] = '\0'; PyObject *o_tmp = 0; - o_tmp = PyObject_GetAttrString(oOpt, "mirSub"); - if(o_tmp == 0) throw strEr_BadOptG; - pOpt->mirSub = ParseSructSRWLOptMirAll(o_tmp, 0, pvBuf, pOpt->mirSubType); - Py_DECREF(o_tmp); - o_tmp = PyObject_GetAttrString(oOpt, "m"); if(o_tmp == 0) throw strEr_BadOptG; if(!PyNumber_Check(o_tmp)) throw strEr_BadOptG; @@ -2525,6 +2701,11 @@ void ParseSructSRWLOptG(SRWLOptG* pOpt, PyObject* oOpt, vector* pvBuf Py_DECREF(o_tmp); } } + + o_tmp = PyObject_GetAttrString(oOpt, "mirSub"); + if(o_tmp == 0) throw strEr_BadOptG; + pOpt->mirSub = ParseSructSRWLOptMirAll(o_tmp, 0, pvBuf, pOpt->mirSubType); + Py_DECREF(o_tmp); } /************************************************************************//** @@ -4204,7 +4385,13 @@ void DeallocOptCntArrays(SRWLOptC* pOptCnt) else if(strcmp(sType, "lens") == 0) delete (SRWLOptL*)(pOptCnt->arOpt[i]); else if(strcmp(sType, "zp") == 0) delete (SRWLOptZP*)(pOptCnt->arOpt[i]); else if(strcmp(sType, "waveguide") == 0) delete (SRWLOptWG*)(pOptCnt->arOpt[i]); - else if(strcmp(sType, "grating") == 0) delete (SRWLOptG*)(pOptCnt->arOpt[i]); + else if(strcmp(sType, "grating") == 0) + { + SRWLOptG *pG = (SRWLOptG*)(pOptCnt->arOpt[i]); + DeallocSRWLOptMir(pG->mirSub, pG->mirSubType); + delete pG; + } + else if(strncmp(sType, "mirror: ", 8) == 0) DeallocSRWLOptMir(pOptCnt->arOpt[i], sType); else if(strcmp(sType, "transmission") == 0) delete (SRWLOptT*)(pOptCnt->arOpt[i]); //{ // SRWLOptT *pT = (SRWLOptT*)(pOptCnt->arOpt[i]); From 4ba7c2cdcd13ef68b02d45ddecd8492104c387ec Mon Sep 17 00:00:00 2001 From: Roman Chernikov Date: Tue, 23 Jun 2026 23:58:06 -0400 Subject: [PATCH 5/6] Helper function for transmission elements --- env/python/srwpy/examples/SRWLIB_Example22.py | 15 ++- env/python/srwpy/uti_mtrl.py | 98 +++++++++++++++++++ 2 files changed, 104 insertions(+), 9 deletions(-) diff --git a/env/python/srwpy/examples/SRWLIB_Example22.py b/env/python/srwpy/examples/SRWLIB_Example22.py index 16414c83..1aa9e936 100644 --- a/env/python/srwpy/examples/SRWLIB_Example22.py +++ b/env/python/srwpy/examples/SRWLIB_Example22.py @@ -18,7 +18,6 @@ from srwpy.uti_plot import * from srwpy.uti_mtrl import * -import math import os from copy import deepcopy @@ -42,6 +41,10 @@ # Aluminum filter filter_thick = 50.e-06 # [m] al_delta, al_atten_len = calc_delta_atten_len('Al', photon_energy) +opFilter = srwl_opt_setup_transm_from_material( + 'Al', filter_thick, photon_energy, + _rx=2.e-03, _ry=2.e-03, _nx=2, _ny=2, _ext_tr=1 +) # B4C mirror; B4C is a formula not registered with a default xraydb density b4c_density = 2.5 # [g/cm^3] @@ -128,14 +131,8 @@ #********************** Optical Elements -# Uniform Aluminum filter. Amplitude transmission is exp(-thickness/(2*L)). -filter_amp = math.exp(-0.5*filter_thick/al_atten_len) -filter_opd = -al_delta*filter_thick -filter_ar_tr = array('d', [filter_amp, filter_opd]*4) -opFilter = SRWLOptT( - _nx=2, _ny=2, _rx=2.e-03, _ry=2.e-03, - _arTr=filter_ar_tr, _extTr=1 -) +# Uniform Aluminum filter. +filter_amp = opFilter.arTr[0] mirror_len = 0.2 mirror_width = 5.e-03 diff --git a/env/python/srwpy/uti_mtrl.py b/env/python/srwpy/uti_mtrl.py index ab1a4d72..ce92a306 100644 --- a/env/python/srwpy/uti_mtrl.py +++ b/env/python/srwpy/uti_mtrl.py @@ -7,6 +7,7 @@ calc_coated_refl_arr calc_multilayer_refl_arr calc_delta_atten_len + srwl_opt_setup_transm_from_material add_mat @@ -22,6 +23,11 @@ import numpy as np +try: + from srwlib import SRWLOptT +except Exception: + from .srwlib import SRWLOptT + xraydb = None xraydb_found = False @@ -478,6 +484,98 @@ def calc_delta_atten_len(_mat, _ph_en, _dens=None): return array('d', delta.ravel()), array('d', atten_len_m.ravel()) +def srwl_opt_setup_transm_from_material( + _mat, + _thick, + _ph_en, + _dens=None, + _rx=1.e-03, + _ry=1.e-03, + _nx=2, + _ny=2, + _x=0, + _y=0, + _ext_tr=1, + _fx=1.e+23, + _fy=1.e+23 +): + """Set up a uniform transmission element from material data. + + :param _mat: material name or chemical formula + :param _thick: material thickness [m] + :param _ph_en: photon energy [eV], or a sequence of photon energies + :param _dens: material density [g/cm^3]; if omitted, use the xraydb material database + :param _rx: horizontal coordinate range [m] + :param _ry: vertical coordinate range [m] + :param _nx: number of points vs horizontal position + :param _ny: number of points vs vertical position + :param _x: horizontal transverse coordinate of center [m] + :param _y: vertical transverse coordinate of center [m] + :param _ext_tr: transmission outside the grid/mesh is zero (0), or same as boundary (1) + :param _fx: estimated focal length in the horizontal plane [m] + :param _fy: estimated focal length in the vertical plane [m] + :return: SRWLOptT by default + """ + delta, atten_len = calc_delta_atten_len(_mat, _ph_en, _dens) + + try: + thick = float(_thick) + if (not np.isfinite(thick)) or (thick < 0): + raise ValueError + except (TypeError, ValueError): + print("Warning: material thickness is invalid. Transmission element is set to vacuum.") + thick = 0.0 + + try: + energy = np.asarray(_ph_en, dtype=float).ravel() + if energy.size < 1: + raise ValueError + except (TypeError, ValueError): + print("Warning: photon energy is invalid. Transmission element mesh energy is set to 0.") + energy = np.asarray([0.0]) + + ne = int(energy.size) + try: + nx = int(_nx) + ny = int(_ny) + if (nx <= 0) or (ny <= 0): + raise ValueError + except (TypeError, ValueError): + print("Warning: transmission mesh dimensions are invalid. A 2 x 2 mesh is used.") + nx = 2 + ny = 2 + + delta_arr = np.asarray(delta, dtype=float).ravel() + atten_len_arr = np.asarray(atten_len, dtype=float).ravel() + if (delta_arr.size != ne) or (atten_len_arr.size != ne): + if (delta_arr.size == 1) and (atten_len_arr.size == 1): + delta_arr = np.full(ne, delta_arr[0]) + atten_len_arr = np.full(ne, atten_len_arr[0]) + else: + print("Warning: material data size is inconsistent with photon energy mesh. Transmission element is set to vacuum.") + delta_arr = np.zeros(ne) + atten_len_arr = np.full(ne, 1.e+23) + + if (not np.all(np.isfinite(delta_arr))) or (not np.all(np.isfinite(atten_len_arr))) or np.any(atten_len_arr <= 0): + print("Warning: material data contain invalid values. Transmission element is set to vacuum.") + delta_arr = np.zeros(ne) + atten_len_arr = np.full(ne, 1.e+23) + + amp = np.exp(-0.5*thick/atten_len_arr) + opd = -delta_arr*thick + ar_tr_one_point = np.empty(2*ne, dtype=float) + ar_tr_one_point[0::2] = amp + ar_tr_one_point[1::2] = opd + ar_tr = array('d', np.tile(ar_tr_one_point, nx*ny)) + + op_t = SRWLOptT( + _nx=nx, _ny=ny, _rx=_rx, _ry=_ry, _arTr=ar_tr, _extTr=_ext_tr, + _Fx=_fx, _Fy=_fy, _x=_x, _y=_y, _ne=ne, + _eStart=float(energy[0]), _eFin=float(energy[-1]) + ) + + return op_t + def add_mat(name, formula, density, categories=None): """Add a material to the user-local xraydb material database.""" if add_material is None: From df81d104913d0eb6f41a32fec6280cd43a0c5d13 Mon Sep 17 00:00:00 2001 From: Roman Chernikov Date: Thu, 25 Jun 2026 09:48:28 -0400 Subject: [PATCH 6/6] Update Pt coating thickness from 200 nm to 20 nm --- env/python/srwpy/examples/SRWLIB_Example22.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/env/python/srwpy/examples/SRWLIB_Example22.py b/env/python/srwpy/examples/SRWLIB_Example22.py index 1aa9e936..4fc9f62e 100644 --- a/env/python/srwpy/examples/SRWLIB_Example22.py +++ b/env/python/srwpy/examples/SRWLIB_Example22.py @@ -1,7 +1,7 @@ ############################################################################# # SRWLIB Example # 22: Material-dependent reflective and refractive optics # Simulating a Gaussian X-ray beam passing through an Aluminum filter, -# a B4C mirror, a 200 nm Pt-coated Si mirror, and an array of diamond CRLs +# a B4C mirror, a 20 nm Pt-coated Si mirror, and an array of diamond CRLs # v 0.01 ############################################################################# @@ -55,8 +55,8 @@ _ang_start=ang_start, _ang_fin=ang_fin, ) -# 200 nm Pt coating on a Si substrate; xraydb expects thickness in Angstroms -pt_coating_thick = 2000. # [Angstrom] +# 20 nm Pt coating on a Si substrate; xraydb expects thickness in Angstroms +pt_coating_thick = 200. # [Angstrom] pt_si_refl = calc_coated_refl_arr( 'Pt', pt_coating_thick, 'Si', _n_ph_en=n_ph_en, _n_ang=n_ang, _n_comp=n_comp,